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INTRODUCTION 


The  objectives  of  this  report  are  to 

(1)  present  algorithms  suitable  for  mini¬ 
computer  application  which  will  per¬ 
form  three  different  general  opera¬ 
tions  used  in  computerized  mapping: 
the  transformation  of  strings  of  con¬ 
tour  data  into  digital  elevation 
models,  the  transformation  of  a  digi¬ 
tal  elevation  model  into  a  set  of 
contour  strings,  and  the  smoothing  of 
a  digital  elevation  model  by  use  of 
various  filtering  and  convolution 
techniques . 

(2)  show  sample  results  of  the  above  op¬ 
erations  applied  to  specific  data 
sets  using  a  range  of  operational 
parameters. 

(3)  report  on  a  comparison  of  the  perfor¬ 
mance  of  two  interpolation  algorithms 
which  can  be  used  to  resample  digital 


elevation  modules 


This  work  was  performed  by  ZYCOR,  Inc. 
under  Contract  DAAK  70-80-C-0248  with  the  U.S.  Army  Engi¬ 
neer  Topographic  Laboratories  ( ETL )  and  the  Defense  Map¬ 
ping  Agency  Hydrographic/Topographic  Center  (DMAHTC). 
The  motivation  for  the  study  was  the  realization  that 
much  of  the  computer  software  now  available  to  DMAHTC  and 
ETL  to  perform  automatic  cartographic  tasks  is  limited  to 
main-frame  computers  and  does  not  make  use  of  the  latest 
developments  in  the  field  of  automated  cartograhy.  If 
the  future  work  of  those  agencies  is  to  be  performed 
using  digital  techniques,  then  superior  algorithms  will 
have  to  be  developed  to  match  the  more  powerful  hardware 
tools  which  are  becoming  available.  With  this  in  mind 
new  algorithms  were  created  to  perform  the  tasks  listed 
under  (1)  above.  With  further  testing  and  development, 
these  algorithms  may  provide  a  software  package  which  can 
handle  most  of  the  different  types  of  data  and  mapping 
problems  to  which  they  are  applicable. 

While  these  algorithms  were  being  develop¬ 
ed  it  was  realized  that  very  little  was  known  about  the 
performance  of  the  available  algorithms  for  carrying  out 
interpolation  operations  on  square  grids  of  two-dimen¬ 
sional  data.  Therefore  part  of  this  study  effort  was  ex¬ 
panded  to  include  a  comparison  of  two  algorithms  which 
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appear  reasonable  candidates  for  handling  digital  terrain 
models  efficiently  and  accurately.  Unlike  the  algorithms 
described  above,  the  grid  resampling  techniques  are  not 
in  themselves  mapping  tools.  Rather,  they  are  algorithms 
which  can  be  used  as  a  part  of  a  complete  grid-to-grid 
resampling  program. 

Two  interim  technical  reports  have  already 
been  delivered  under  this  contract.  The  first  discussed 
the  mathematical  development  and  preliminary  design  con¬ 
siderations.  The  second  (Reference  1)  discussed  program 
inputs,  outputs,  and  control  parameters,  and  also  pre¬ 
sented  some  sample  results.  Material  in  '.hose  reports 
relevant  to  the  development  of  the  algorithms  is  included 
in  the  final  report. 

Section  2  of  this  report  describes  the 
test  data  made  available  by  the  Engineer  Topographic  Lab¬ 
oratories  ( ETL ) .  This  data  was  originally  in  the  form  of 
contour  strings  covering  a  7.5  minute  rectangle  in  south¬ 
ern  Arizona.  It  was  broken  .nto  strings  corresponding  to 
smaller  areas  and  used  in  several  different  ways  to  pro¬ 
vide  data  for  testing  the  algorithms  discussed  in  this 
report . 


Section  3 

discusses  the  algorithm 

to 

per- 

form 

the  transformation 

of  contour 

strings  and 

sim 

ilar 

data 

such  as 

drain  and 

ridge 

lines 

and  lake  boundaries 

into 

dig ital 

terrain  models. 

The 

development 

of 

the 
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algorithm  is  described.  An  overview  of  its  operation  is 
given.  Results  of  testing  this  algorithm  using  the  data 
described  in  Section  2  are  provided. 

Section  4  discusses  the  algorithm  to  per¬ 
form  the  transformation  of  digital  terrain  models  into 
contour  strings.  This  algorithm  is  essentially  the 
inverse  of  the  contour- to-grid  algorithm.  The  develop¬ 
ment  of  the  algorithm  is  described.  An  overview  of  its 
operation  is  given.  Again,  results  of  testing  this  algo¬ 
rithm  using  the  data  described  in  Section  2  are  provided. 

Section  5  discusses  a  variety  of  algo¬ 
rithms  for  smoothing  terrain  elevation  models.  Use  of 
these  algorithms  may  be  an  independent  step  or  a  post 
processing  step  which  is  part  of  the  grid-to-contour  or 
contour- to-grid  algorithms.  Techniques  considered  in¬ 
clude  one  step  least  squares  operators  and  recursive  con¬ 
volution  operators  with  a  number  of  different  con¬ 
straints.  Testing  was  carried  out  and  conclusions  and 
recommendations  are  provided. 

Section  6  is  a  comparison  of  the  Akima  and 
Jancaitis  algorithms  for  two-dimensional  interpolation. 
They  are  compared  based  on  statistical  performance,  visu¬ 
al  interpretation,  timing,  and  frequency  characteris¬ 
tics.  Recommendations  are  made  on  choosing  between  the 
two  algorithms. 
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Finally,  Section  7  contains  a  general  re¬ 
view,  conclusions,  and  recommendations,  including  sum¬ 
maries  of  those  given  in  prior  sections.  A  discussion  of 


future  work  required  in  this  area  is  included. 
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2. 


TEST  DATA 


DMA/ETL  provided  2YCOR  with  digitized  con¬ 
tour  data  for  a  7.5  minute  x  7.5  minute  quadrangle  in 
southern  Arizona  to  be  used  as  the  basis  for  test  data 
for  this  report.  Also  provided  were  digitized  drain  and 
ridge  lines  and  digitized  lake  boundaries  for  the  same 
area.  All  data  were  digitized  at  a  resolution  of  approx¬ 
imately  16  meters. 

The  supplied  data  had  problems  which  re¬ 
quired  attention.  In  areas  of  large  gradient  the  auto¬ 
matic  digitization  process  utilized  to  create  the  contour 
strings  from  an  existing  contour  map  had  difficulty  in 
following  individual  contours.  The  result  was  an  occa¬ 
sional  pattern  of  digitized  contours  which  crossed  each 
other  or  terminated  suddenly  and  arbitrarily  .  The  pres¬ 
ence  of  erroneous  data  made  it  difficult  to  find  suitable 
test  areas,  especially  ones  with  significant  terrain 
variation.  Careful  picking  of  test  areas  surmounted  this 
problem.  The  drain  and  ridge  data  were  extremely  noisy. 
Editing  of  the  data  was  required  both  to  smooth  it  and  to 
make  it  agree  with  the  contour  data. 

The  data  were  partitioned  by  ZYCOR  into 
smaller  test  areas  which  provided  inputs  to  the  various 
programs  described  in  this  report.  The  raw  contours, 


fa 


drains  and  ridge  lines,  and  lake  boundaries  were  used  as 
inputs  to  the  contour-to-gr id  algorithm  discussed  in  Sec¬ 
tion  3.  The  digitized  terrain  models  produced  were  used 
as  inputs  to  the  gr id- to-contour  algorithms  discussed  in 
Section  4,  and  to  the  smoothing  and  generalization  algo¬ 
rithms  discussed  in  Section  5.  Figures  2.1  and  2.2  show 
the  digitized  areas  which  provided  the  test  sets  utilized 
in  Sections  3  and  5.  The  labels  along  the  sides  corres¬ 
pond  to  the  digitizing  coordinates. 

Three  small  local  areas  were  chosen  to  be 
used  in  testing  the  gr id- to-con tour  algorithm  discussed 
in  Section  4.0.  Other  larger  areas  were  utilized  as  data 
for  the  resampling  study  as  discussed  in  Section  6.  The 
test  areas  utilized  are  discussed  in  that  section. 


7 


7<*0  72Q  700  680  660 


09*  09*  DOS  OZS  0*S 


CONTOUR-TO-GRID  INTERPOLATION 


3.1  OVERVIEW.  Contour- to-Gr id  (CTOG)  is  the  name 
of  a  new  computer  program  which  transforms  vectorized 
contours  and  related  geomorphic  data  into  digital  terrain 
models  ( DTM )  .  It  is  expected  that  CTOG  will  become  an 
important  step  in  DMA  operations  to  produce  a  digital 
data  base  from  its  extensive  library  of  maps  and  charts. 
With  some  added  work,  the  CTOG  algorithm  could  become  a 
powerful  grid  editing  tool.  In  that  capacity  it  would  be 
used  to  synthesize  elevation  values  from  contour  maps  in 
areas  where  photogrammetr ic  techniques  were  not  satisfac¬ 
tory. 

The  CTOG  algorithm  is  made  up  of  three 
processing  steps.  First,  initial  estimates  of  the  grid 
node  values  are  obtained  by  interpolation  from  the  input 
data.  Second,  those  grid  nodes  which  are  "close"  to  the 
input  data  are  adjusted  to  fit  the  data  and  fixed  so  that 
they  will  not  be  changed  by  the  last  processing  step. 
Finally,  the  initialized,  adjusted,  and  fixed  grid  is 
filtered  by  a  spatial  convolution  operator  to  produce  a 
smoothly  varying  surface  in  regions  not  closely  tied  to 
the  input  data.  An  overview  of  these  thee  steps  is  givn 
in  Section  3.2  to  3.4  below.  Certain  technical  details 
are  included  in  the  appendices.  A  description  of  test 
procedures  and  test  results  is  given  in  Section  3.5. 
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3.2  IN  IT I AL I Z AT ION .  The  initialization  techniques 
employ  contours,  ridge  and  streai.  lines,  lake  boundaries 
and  point  elevations  to  derive  a  first  approximation  of 
element  values  at  the  nodes  of  the  grid.  The  initializa¬ 
tion  algorithm  attempts  to  make  use  of  the  highly  struc¬ 
tured  nature  of  the  input  data  in  a  manner  similar  to 
that  used  by  an  individual  creating  the  grid  by  hand  from 
a  contour  map  would  make.  To  do  this  the  algorithm1  uses 
the  contour  lines  which  the  connected  digitized  points 
represent  rather  than  just  the  points  themselves.  This 
contrasts  strongly  with  conventional  gridding  algorithms 
which  usually  treat  the  digitized  contour  points  as 
random  elevation  data. 

Figure  3.1  illustrates  the  method  for 
selecting  the  contour  lines  which  are  used  to  compute  an 
initial  estimate  for  a  grid  node  elevation.  In  this  par¬ 
ticular  example,  four  data  points  are  selected  to  be  used 
to  estimate  the  elevation  at  the  indicated  node.  The 
four  points  are  found  by  searching  two  pairs  of  direc¬ 
tions  (up-down,  and  left-right  in  this  case)  until  con¬ 
tour  lines  are  found.  A  point  is  selected  at  the  first 
intersection  of  a  contour  line  with  each  of  the  straight 
lines  that  emanate  from  the  grid  node. 


[3  zycoR  INC 


FIGURE  3.1 


Selection  of  Points  Used  to  Estimate,  Elevation  at 
Grid  Node.  (The  number  of  points  and  the  search 
directions  are  not  necessarily  limited  to  those  shown.) 
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By  looking  at  Figure  3.1,  it  is  easy  to 
see  that  interpolation  of  the  node  elevation  by  using  the 
pair  of  contour  intersections  provided  by  the  horizontal 
search  lines  is  the  method  which  would  be  used  by  an  in¬ 
dividual  performing  this  task  by  hand.  This  algorithm 
obviously  cannot  duplicate  the  performance  of  a  human 
being  performing  this  task,  but  it  comes  close  by  using 
many  search  directions,  up  to  8  including  the  4  in  Figure 
3.1  and  4  more  rotated  it  a  45°  angle  to  those  shown;  and 
by  weighting  each  pair  of  contour  intersections  along 
opposite  search  directions  by  distance  or  slope  dependent 
weighting  function.  The  final  estimate  of  each  grid  node 
elevation  is  the  weighted  sum  of  each  estimate  obtained 
by  using  a  different  pair  of  search  directions. 

The  weighting  functions  can  be  chosen 
either  to  emphasize  the  search  directions  which  corres¬ 
pond  to  short  distances,  or  to  emphasize  large  differ¬ 
ences  between  intersected  contours.  In  the  example 
shown,  both  choices  of  weighting  functions  would  give  the 
highest  influence  to  the  horizontal  search  direction  pair 
which  is  the  correct  choice.  Weighting  function  choices 
are  discussed  in  Appendix  E.  Interpolation  between  con¬ 


tour  intercepts  may  be  either  linear  or  quadratic. 


In 


the  latter  case,  it  is  necessary  to  find  two  contour 
intersections  in  each  search  direction  rather  than  one. 


Interpolation  schemes  are  discussed  in  Appendix  A 

Figure  3.2  shows  the  more  complicated  case 
using  all  8  search  directions.  This  figure  indicates  the 
principle  complexities  encountered  in  developing  this 
algorithm;  efficiently  finding  and  using  the  contour  in¬ 
terceptions  with  search  lines.  This  is  done  by  first 
searching  along  each  contour  for  the  search  line  inter¬ 
cepts  and  then  sorting  the  intercepts  for  each  search 
line  so  that  they  may  be  easily  found  starting  at  any 
grid  node. 

In  this  init ialization  step,  digitized 
drain  and  ridge  lines  as  well  as  lake  boundaries  may  be 
utilized  in  a  manner  similar  to  digitized  contours. 
Unlike  contours  and  lake  boundaries  which  have  a  single 
elevation  associated  with  every  point,  a  drain  or  ridge 
line  has  a  different  elevation  at  each  point.  Thus  the 
elevation  at  the  intersection  of  a  search  line  with  one 
of  these  auxiliary  inputs  must  be  computed  by  linear  in¬ 
terpolation  bvetween  digitized  points. 

3.3  ADJUSTMENT  AND  FIXING.  After  initialization 
the  grid  adjustment  step  may  be  used  to  force  the  surface 
to  pass  through  given  control  points.  The  control  points 


Figure  3.2.  Illustrations  of  Eight  "Search  Lines"  Used 
to  Determine  Interpolation  Points  for  Each 
Grid  Node.  (The  horizontal  and  vertical 
search  lines  are  the  same  as  the  grid  lines 
of  the  target  matrix.  The  diagonal  search 
lines  are  dashed  to  distinguish  them  from 
grid  lines.  Grid  nodes  are  shown  as  small 
dots,  and  interpolation  points  (where  search 
lines  intersect  contours)  are  shown  as  large 
dots . ) 
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include  the  intersections  of  grid  lines  with  contours, 
ridge  and  drain  lines  and  lake  boundaries.  They  may  also 
include  spot  elevations.  The  grid  nodes  within  a  user 
chosen  distance  of  the  control  points  are  fixed  so  that 
subsequent  processing  will  not  change  their  values. 
Also,  all  grid  nodes  interior  to  lake  boundaries  are 
fixed  at  this  time. 

The  adjustment  step  involves  fitting  a 
surface  of  the  form  axy  +  bx  +  cy  +  d  to  the  four  grid 
nodes  surrounding  a  control  point  and  then  raising  or 
lowering  the  surface  until  it  passes  through  the  control 
point.  Evaluation  of  the  surface  at  the  grid  nodes  pro¬ 
vides  the  new  fixed  values. 

This  step  is  primarily  important  in  cases 
in  which  it  is  desired  to  validate  the  algorithm  by  re¬ 
contouring  the  output  grid  at  levels  matching  those  of 
the  input  contour  strings.  The  adjustment  and  fixing 
guarantees  that  the  contours  corresponding  to  input  con¬ 
tour  strings  will  not  be  significantly  shifted  anywhere 
in  the  entire  map. 

Discuosions  with  DMA  personnel  have  indi¬ 
cated  that  this  is  not  a  normal  goal.  Thus  grid  adjust¬ 
ment  may  be  a  step  the  user  will  often  wish  to  skip. 
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3.4  FILTERING.  The  final  processing  step  involves 

applying  a  convolution  filter  to  the  initial  grid  to  gen¬ 


erate  a  smoothly  varying  gridded  surface.  This  will 
eliminate  much  of  the  roughness  of  the  grid  which  is 
caused  by  using  local  operations  to  perform  the  initiali¬ 
zation  process.  All  nodes  of  the  grid  are  allowed  to 
change  values  under  this  filtering  process,  except  those 
fixed  because  they  were  close  to  the  input  contours  or 
other  data  structures. 

The  convolution  operator  used  in  this  step 
is  either  a  thirteen  point  bi-harmonic  operator  or  a  5 
point  Laplacian  operator.  Filtering  continues  for  a 
specified  number  of  cycles  or  until  the  fractional  change 
in  grid  values  between  cycles  is  below  a  set  value. 
These  spatial  filtering  algorithms  are  the  same  ones  used 
in  the  grid  smoothing  operations  discussed  in  Section  5 
of  this  report.  More  information  on  their  characteris¬ 
tics  and  methods  of  terminating  their  operations  are 
available  in  that  Section. 

3.5  CONTOUR  TO  GRID  TESTING.  Testing  of  this  algo¬ 
rithm  was  carried  out  using  the  DMA  supplied  data  as 
described  in  Section  2.  Evaluation  of  results  must  be 
subjective  since  ground  truth  data  does  not  fxist.  Thus, 
in  this  section  as  in  those  which  follow,  visual  results 


will  be  emphasized.  Two  different  types  of  output  are 
presented:  first,  contour  maps  created  from  grids  pro¬ 

duced  by  CTOG  which  demonstrate  an  ability  to  replicate 
the  original  contour  inputs  and  to  handle  auxiliary  in¬ 
puts  such  as  drains,  ridges  and  lake  boundaries;  and 
second,  very  small  posted  areas  in  which  the  grid  node 
values  calculated  can  be  judged  against  the  raw  data. 

Figures  3.3  and  3.4  are  contour  maps  pro¬ 
duced  from  CTOG  runs  using  the  two  test  areas  described 
in  Sections  2.0  and  shown  in  Figures  2.1  and  2.2.  One 
pass  of  the  bi-harmonic  filter  was  applied  to  the  65x65 
output  grids  in  these  runs.  The  effects  of  additional 
filtering  on  these  figures  will  be  discussed  in  Section 
5.0  of  this  report. 

Nodes  fixed  by  the  adjustment  step  of  CTOG 
are  indicated  in  Figures  3.3  and  3.4  by  small  x's.  The 
large  area  of  fixed  nodes  in  the  lower  left  corner  of 
Figure  3.3  corresponds  to  the  interior  of  a  lake.  All 
grid  nodes  within  the  lake  boundary  are  fixed  to  the 
level  on  the  lake  boundary. 

For  both  maps  the  match  between  the  con¬ 
toured  CTOG  output  and  corresponding  input  data  is  good. 
However,  the  closures  in  the  upper  part  of  3.3  did  not 
match  the  original.  Part  of  the  difference  is  attribut¬ 
able  to  the  contouring  logic  used  to  create  the  maps  in 
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Figure  3.4.  CTOG  Output  for  Test  Area  B. 
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this  section.  Post  processing  using  one  of  the  filtering 
options  would  also  improve  these  features  (see  Section 
5). 

Figures  3.5  and  3.6  show  small  subsections 
of  3.3  and  3.4.  The  grid  node  positions  are  shown  by  the 
+'s  and  their  elevation  values  are  posted.  The  lines  re- 
present  the  input  contour  strings  and  drain  lines.  Fixed 
nodes  are  indicated  by  star  shaped  figures. 

Figure  3.5  corresponds  to  part  of  the  low¬ 
er  left  corner  of  3.3.  In  it  can  be  seen  the  lake  bound¬ 
ary  surrounding  the  dense  pattern  of  fixed  nodes  all  with 
the  same  posted  level  of  4860.  A  few  nodes  outside  of 
the  lake  are  also  fixed  during  the  adjustment  step. 

Figure  3.6  corresponds  to  an  area  in  the 
right  center  of  3.4.  It  provides  a  detailed  view  of  the 
way  fixed  nodes  are  clustered  along  input  data.  For 
these  figures  all  nodes  closer  than  1/2  a  cell  diagonal 
to  a  control  point  are  fixed.  As  defined  in  Section  3.3 
control  points  include  all  points  used  to  define  the 
input  contours  plus  the  crossing  of  grid  rows  and  columns 
by  contour  strings. 
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Figure  3.5.  Detail  of  CTOG  Output, 
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Figure  3.6.  Detail  of  CTQG  Output. 
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4. 


GRID-TO-CONTOUR 


4.1  OVERVIEW.  Grid-to-contour  (GTOC)  is  the  name 
used  to  identify  a  new  DMA  computer  program  which  inter¬ 
polates  elevation  contours  from  digital  terrain  models 
( DTM ) .  That  is,  the  major  input  for  the  program  is  a  DTM 
while  the  output  is  a  set  of  vectorized  contours  that  can 
be  plotted  on  several  of  the  graphic  systems  at  DMA. 

An  immediate  application  of  GTOC  in  DMA  is 
to  analyze  and  verify  DTM 1 s  produced  by  the  contour-to- 
grid  (CTOG)  program  and  other  sources  in  DMA.  Its  long 
term  applications  will  probably  be  in  map  production  us¬ 
ing  as  source  data  the  extensive  library  of  DTM's  being 
assembled  by  DMA. 

DMA  has  developed  similar  contour  genera¬ 
tion  programs.  They  are  confined  to  the  mainframe  com¬ 
puting  system  and  are  used  primarily  for  data  verifica¬ 
tion  and  evaluation.  The  new  technology  and  design  of 
GTOC  will  enable  DMA  to  use  the  program  on  many  of  the 
mini-computer  based  systems  currently  online  or  planned 
for  DMA.  Also,  the  quality  of  output  suggests  the  pro¬ 
gram  will  be  useful  in  future  map  production  plans  for 
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DMA. 
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A  description  of  the  logic  and  major  algo¬ 
rithms  used  in  GTOC  are  provided  in  the  following  sec¬ 
tions.  Several  test  outputs  based  on  DTM's  produced  by 
CTOG  are  provided  as  examples  of  the  program  outputs. 
Further  testing  and  evaluation  of  GTOC  is  being  performed 
by  ETL  and  DMA  personnel. 

4.2  BACKGROUND.  The  general  contouring  problem  is 
to  establish  all  points  (x,y)  which  satisfy  the  equation 

E(x,y)  -C  *  0 

where  E  represents  the  elevation  as  it  varies  with  x  and 
y  (easting  and  northing  coordinates)  and  C  is  the  contour 
val ue . 

Cartographers  "solve’'  this  equation  using 
a  variety  of  tools,  logic,  and  experience.  Their  tools 
range  from  simple  geometric  dividers  to  complex  analog 
photogrammetry  systems.  The  data  they  use  may  be  a  set 
of  widely  spaced  spot  elevation  data  or  continuous  photo¬ 
graphic  images  of  the  topography. 

Various  digital  systems  used  by  the  gov¬ 
ernment  require  elevation  data  in  a  form  that  is  easy  to 
handle  while  maintaining  high  accuracy,  nearly  equivalent 
to  the  analog  images  of  E(x,y).  This  requirement  is 
satisfied  by  digital  terrain  models  which  are  matrices  of 
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elevation  values.  The  values  correspond  to  points  along 
a  series  of  equally  spaced  profiles  over  the  x-y  plane. 
Points  along  the  profiles  are  also  equally  spaced. 

Although  cartographers  could  apply  their 
classical  contouring  methods  to  the  tabulated  elevations 
in  a  DTM,  this  is  unrealistic  because  of  the  volume  of 
data  and  the  rate  at  which  DTM's  can  be  produced.  Hence 
computer  contouring  is  a  reasonable  alternative. 

In  its  simpliest  form,  computer  contouring 
is  a  process  of  connecting  dots.  For  example,  imagine 
that  the  DTM  representing  E(x,y)  is  a  grid  of  elevations 
where  vertical  grid  lines  correspond  to  adjacent  profiles 
of  E(x,y)  and  horizontal  lines  connect  values  on  adjacent 
profiles.  If  a  contour  level  lies  between  two  DTM 
values,  then  the  crossing  point  is  on  the  grid  line  con¬ 
necting  the  values.  Other  dots  can  be  positioned  on  grid 
lines  where  this  contour  crosses.  Finally,  the  contour 
curve  can  be  drawn  by  connecting  the  dots  with  short  line 
segments  or  curves.  Of  course,  there  is  an  order  that 
must  be  followed  in  connecting  the  dots.  Figure  4.1 
illustrates  these  simple  steps  for  a  small  portion  of  a 
DTM. 
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Contours  Drawn  by  connecting  Contouis  Drawn  by  connecting 

the  Points  with  Lane  Segments  the  Kants  with  Smooth  Curves 


Figure  4.1.  Ccniputei  Contoui  ing  using  Lane  Segments  or  Smooth  Cuives 


The  accuracy  and  acceptability  of  contours 
produced  by  this  simple  process  depend  on  several  fac¬ 
tors.  These  are 

1.  accuracy  of  the  DTM  values, 

2.  spacing  between  DTM  values, 

3.  texture  of  the  area  covered  by  the 
DTM, 

4.  procedure  for  computing  dot  locations, 
and 

5.  procedures  for  connecting  the  dots. 
Only  Factors  4  and  5  are  within  the  control  of  the  con¬ 
touring  logic  described  above.  It  assumes  that  the  DTM 
data  are  accurate  and  properly  spaced  for  the  topography 
it  represents.  Indeed,  Factors  4  and  5  are  the  reasons 
why  computer  contouring  is  not  as  simple  as  suggested  in 
the  discussion  above.  There  are  a  variety  of  techniques 
for  each.  Several  were  investigated  and  the  most  effec¬ 
tive  combination  implemented  for  DMA.  These  are  discuss¬ 
ed  in  the  next  sections. 

4.3  CONTOURING  METHODS.  Two  methods  for  implement¬ 

ing  contouring  logic  were  tested.  One  was  selected  as 
the  more  desirable  and  installed  at  DMA.  Both  are  de¬ 
scribed  here. 
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The  first  method,  called  the  "cell-prior¬ 
ity,"  was  initially  proposed  to  DMA,  implemented,  and 
tested.  Because  of  problems  discussed  below,  it  was  re¬ 
placed  by  the  "curve-priority"  method  which  is  discussed 
in  Section  4.3.2. 


4.3.1  Cell-Priority  Contouring.  Cell  priority  con¬ 
touring  works  as  follows:  the  DTM  grid  lines  define  a 
collection  of  grid  cells,  each  bounded  by  adjacent 
horizontal  and  vertical  grid  lines.  There  are  elevation 
values  at  the  corner  of  each  cell.  Processing  on  the 
cells  starts  at  the  upper  left  corner  of  the  DTM  and 
steps  down  and  across  the  area.  The  first  cell's  corner 
values  are  selected  and  compared  against  all  contour 
values  to  be  drawn.  If  a  contour  intersects  the  sides  of 
the  cell,  the  dot  locations  are  computed  and  stored  with 
the  cell.  If  several  curves  pass  through  the  cell,  then 
there  will  be  a  set  of  dots  for  each  curve,  each  set 
stored  with  the  cell.  Once  all  curves  through  a  cell  are 
stored,  the  next  cell  is  retrieved  and  processed.  This 
process  continues  until  all  cells  in  the  DTM  have  been 
processed. 

When  the  cell  size  at  map  scale  is  very 
small,  connecting  the  two  dots  that  define  a  curve  across 
the  cell  might  be  acceptable.  However,  for  lar^r  cells. 
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unrealistic  and  generally  unacceptable  polygonal  shaping 
results.  Consequently,  an  algorithm  was  designed  to 
compute  intermediate  dots  along  the  contour  wiinin  the 
grid  cells.  By  making  the  dots  close  enough,  lines  can 
be  used  to  smoothly  connect  the  dots.  The  method  for 


computing  intermediate  dots 

is  the  same 

for 

cell-  or 

curve-priority 

contouring ; 

hence,  it 

is 

covered 

separately  in 

Section  4.3.4. 

It  suffices 

to 

say  that 

when  it  is  employed,  the  number  of  points  representing 
the  curve  can  increase  from  a  minimum  of  2  to  30  or  more 
depending  on  the  curve's  complexity  across  the  cell.  All 
of  these  points  must  be  stored  with  the  cell  or  displayed 
immediately. 

Cell-priority  contouring  is  very  efficient 
in  some  respects.  Once  a  cell  has  been  processed,  the 
computer  memory  occupied  by  its  corner  values  and  contour 
dots  can  be  released.  No  processing  on  that  cell  will 
occur  again.  Consequently,  small  portions  of  very  large 
DTM's  can  be  loaded,  processed,  and  forgotten.  If  neces¬ 
sary,  dot  locations  on  each  curve  through  a  cell  can  be 
stored  on  discs.  Because  of  this,  the  method  is  very  at¬ 
tractive  for  use  on  mini-computers  with  limited  memory 
and  cheap  disc  storage.  The  method  is  even  more  attrac¬ 
tive  if  the  curves  can  be  displayed  immediately  and  never 
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stored.  On  a  CRT  graphics  device,  the  output  appears  to 
flow  across  the  screen  from  le f t- to-r ig ht ,  filling  the 
screen  with  narrow  vertical  strips  of  the  finished  con¬ 
tour  plot.  Since  there  are  a  variety  of  mini-computers 
and  raster  type  display  devices  at  DMA,  it  was  reasonable 
to  pursue  this  method  for  contouring. 

Problems  develop  when  the  curves  must  be 
drawn  by  an  incremental  plotter,  annotated,  edited,  or 
adjusted.  All  are  important  operations  at  DMA.  Each 
problem  is  a  consequence  of  the  cell-by-cell  processing. 

Since  points  on  curves  are  not  stored  in 
order  along  the  curves  (they  are  stored  by  cells),  it  was 
necessary  to  link  the  segments  together  tor  plotting  or 
for  the  other  operations.  Curve  linking  was  accomplished 
by  processing  each  contour  one-by-one.  When  an  unpro¬ 
cessed  curve  was  located  in  a  cell,  the  points  were  re¬ 
trieved  to  start  a  curve  string.  When  a  curve  passed  to 
an  adjacent  cell,  the  new  cell's  data  were  loaded,  then 
the  points  on  the  curve  were  extracted  and  added  to  the 
string.  This  continued  until  all  cells  containing  points 
on  the  curve  were  loaded  and  added  to  the  string.  The 
resulting  "continuous"  string  of  points  along  the  curve 
were  available  to  plot,  annotate  with  elevation  values, 
or  for  other  operations. 


Although  curve  linking  appears  to  be 
fairly  simple,  the  implementation  within  mini-computer 
constraints  was  difficult  and  not  very  satisfactory. 
Data  volume  became  a  major  problem.  When  curves  are 
defined  by  2  points  per  cell  and  there  is  only  one  curve 
per  cell,  the  amount  of  curve  data  is  usually  smaller 
than  the  original  DTM  data.  However,  when  intermediate 
points  through  cells  are  generated  and  when  there  are 
multiple  curves  in  a  cell  (two  cases  that  must  be  ex¬ 
pected),  the  curve  data  volume  can  rapidly  exceed  the 
storage  requirements  for  the  original  DTM. 

Since  curve  complexity  dictates  the  amount 
of  data  along  a  curve,  it  is  difficult  to  predict  a  pri¬ 
ori  how  much  space  will  be  required  to  load  data  for  each 
curve.  Consequently,  a  variety  of  data  partitioning  and 
blocking  schemes  were  tested.  Since  the  amount  of  com¬ 
puter  memory  was  fairly  restrictive,  none  of  the  schemes 
were  very  satisfactory.  Although  linking  was  achieved, 
it  was  not  finali?ed  because  of  the  operator's  expertise 
required  to  set  control  parameters  and  its  extensive  use 
of  computer  time  and  resources. 

Since  DMA  wanted  the  contouring  program  to 
operate  on  a  variety  of  systems,  mostly  mini-computer 
based,  and  since  the  linking  problem  appeared  to  be  ex¬ 
tremely  time  consuming,  the  cell-priority  method  was  set 
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aside.  The  curve-priority  method  which  partially  solves 
those  problems  was  developed  instead. 

4.3.2  Curve-Priority  Contouring.  This  method  re¬ 

verses  the  processing  priority  on  curves  and  cells. 
Whereas  the  former  method  started  with  a  cell  and  then 
examined  the  cell's  contents  for  all  the  contours,  this 
method  starts  with  a  contour  and  examines  all  the  cells. 
The  processing  starts  with  a  contour  level  C,  initially 
the  first  level  of  interest.  Then,  a  systematic  search 
of  the  cells  is  started  to  locate  a  cell  containing  the 
curve.  The  search  commences  at  the  upper  left  corner  of 
the  matrix  and  scans  down  and  across.  When  a  curve  at 
level  C  is  detected,  points  on  the  curve  through  the  cell 
are  computed  as  described  in  Section  4.3.4.  Since  the 
initial  cell  shares  an  edge  with  a  cell  which  also  con¬ 
tains  the  curve,  it  is  possible  to  move  left,  right,  up, 
or  down  in  a  predictable  manner  from  one  cell  to  the  next 
along  the  path  of  the  curve.  each  time  a  new  cell's  ele¬ 
vation  values  are  retrieved  from  the  matrix,  points  along 
the  curve  are  computed  and  added  to  a  string  representing 
the  contour. 

Tracing  the  curve  from  cell-to-cell  stops 
when  the  last  cell  processed  is  on  the  edge  of  the  matrix 


33 


(contour  exits  on  a  side  of  the  map)  or  when  the  next  new 
cell  is  the  same  as  the  initial  ceil  (the  curve  is  closed 


within 

the 

map)  . 

After 

the 

curve- tracing 

stops , 

the 

next 

curve 

(  if 

any)  at 

level 

C  is  sought. 

The 

searching 

starts 

in 

the  cell 

just 

beyond  the  one 

where 

the 

last 

curve 

was 

initially 

detected.  Searching 

continues 

until 

another  curve  is  located  which  means  that  the  curve-fol¬ 
lowing  logic  is  evoked  again  or  until  the  last  grid  cell 
in  the  matrix  is  analyzed.  The  entire  process  repeats 
for  the  next  contour  value  or  terminates  if  there  are  no 
more  contours. 

To  prevent  the  detection  of  a  single  curve 
more  than  once,  a  logic  matrix  is  constructed.  The  ma¬ 
trix  contains  two  "YES-NO"  flags  per  grid  cell  represent¬ 
ing  the  top  and  left  sides  of  the  cell  area.  Initially, 
for  each  contour  level,  all  flags  are  set  to  "NO".  As  a 
contour  is  traced  through  a  ceil  and  found  to  cross  the 
top  side,  the  top  flag  is  changed  to  "YES".  Similarly 
the  left  side  flag  is  changed  to  "YES"  if  it  crosses  the 
left  side  of  the  cell.  Therefore,  all  cell  areas  where 
the  curve  intersects  the  top  or  left  edges  will  be  marked 
as  having  been  processed.  Later  when  such  a  cell  is 


checked  by  the  curve  detection  algorithm,  it  will  be 
ignored  so  that  the  curve  will  not  be  retraced  a  second, 
third,  etc.  times.  A  new  curve  will  be  started  only  if 
the  top  or  left  side  flags  are  "NO".  Even  if  a  contour 
is  traced  through  both  the  left  side  and  top  of  a  cell, 
this  logic  does  not  prevent  a  second  contour  of  the  same 
level  from  crossing  the  cell  since  it  can  be  picked  up  in 
some  other  cell  and  then  eventually  traced  through  this 
one . 

As  the  contour  strings  are  generated,  they 
are  transferred  directly  to  disc  storage  for  later  pro¬ 
cessing  (display,  editing,  etc.).  The  strings  are  auto¬ 
matically  ordered  by  virtue  of  the  tracing  process. 

Computer  memory  and  data  handling  problems 
are  easier  to  handle  with  the  curve-priority  method. 
Since  the  only  data  in  memory  is  the  DTM,  it  can  easily 
be  partitioned  into  a  number  of  sub-matrices,  each  one  of 
which  will  fit  in  memory.  Contours  within  each  sub-ma¬ 
trix  are  generated  and  output  to  the  disc.  When  the  pro¬ 
gram  is  operated  on  16  bit  mini-computers,  there  are  ap¬ 
proximately  32,000  locations  available  to  store  elevation 
values.  For  a  901x901  DTM,  25  sub-matrices  are  produced, 
5  across  and  5  down  the  area  covered  by  the  DTM.  Con¬ 
tours  are  generated  in  order,  top  left  to  bottom  right. 
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Consequently,  when  plotted,  they  appear  to  be  generated 
in  5  vertical  strips  across  the  area. 


4.3.3  Curve  Drawing  and  Annotation.  Contours  gener¬ 
ated  by  either  of  the  two  methods  are  represented  by 
strings  of  (x,y)  points  along  the  curves.  They  are 
stored  on  discs  for  later  drawing  or  input  to  other  steps 
such  as  editing,  generalization,  etc.  The  only  operation 
provided  for  under  this  effort  was  contour  drawing. 

The  contour  drawing  algorithm  is  fairly 
limited  by  DMA  standards.  It  was  implemented  only  to 
display  and  label  the  curves.  More  specialized  annota¬ 
tion  of  the  contours  is  planned  in  future  work. 

Curves  are  drawn  in  the  order  they  are 
output  by  the  contour  generation  program.  Parameters 
provided  by  the  operator  control  spacing  of  labels  along 
curves  and  label  size.  Additional  parameters  are  pro¬ 
vided  to  select  different  size  drawing  pens  or  line 
widths,  provided  the  plotting  device  is  properly 
equipped. 


4.3.4  Interpolating  Contour  Points.  When  the  area  of 
a  grid  cell  is  very  small  at  map  scale,  contours  can 
probably  be  drawn  satisfactorily  by  connecting  the  points 
on  the  cell  sides  with  line  segments.  However,  when  the 
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cells  are  fairly  large,  it  will  be  necessary  to  interpo- 
late  intermediate  points  through  the  cell  to  produce 
smooth  and  realistic  contours.  The  techniques  used  to 
locate  points  along  a  contour  within  a  cell  are  described 
briefly  here.  More  details  are  provided  in  Appendix  B. 

To  accommodate  both  requirements,  i.e., 
defining  curves  by  edge  intersections  (2  points)  alone  or 
by  multiple  points  across  each  cell,  two  levels  of  com¬ 
plexity  are  required.  Ideally  the  simpliest  requirement 
would  be  satisfied  by  Linear  Fitting  in  which  the  contour 
intersections  with  the  cell  edges  are  joined  by  a  single 
straight  line.  However,  since  some  cells  have  more  than 
two  crossings  of  a  contour  level  with  cell  edges,  a 
further  refinement  is  required  to  prevent  ambiguities. 
Thus  for  this  simpliest  case  a  pattern  of  four  triangles 
is  created,  each  triangle  defined  by  a  cell  side  and  the 
center  of  the  cell.  The  interpolation  is  in  fact  linear 
over  these  triangles. 

Non-Linear  Fitting  yields  two  points  where 
the  curve  intersects  the  cell  edges  plus  one  or  more 
intermediate  points  across  the  cell.  When  these  are 
closely  spaced  and  connected  by  straight  line  segments, 
the  results  appear  to  be  a  smooth  curve  through  the 
cell. 
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Non-Linear  Fitting  is  much  more  complex 
and  time-consuming  since  a  local  elevation  model  must  be 
constructed  across  the  grid  cell.  This  model  is  formu¬ 
lated  from  16  elevation  values  which  are  at  the  corners 
of  a  4x4  sub-matrix  containing  the  grid  cell  at  its  cen¬ 
ter.  Figure  B.3  in  Appendix  B  illustrates  the  extraction 
of  the  4x4  matrix  from  the  DTM . 

Two  algorithms,  called  the  Iterative  and 
Stepping  Methods,  were  developed  to  trace  contours  across 
the  local  elevation  models.  The  Iterative  Method  is  used 
most  frequently  while  the  Stepping  Method  is  used  only 
when  the  Iterative  Method  fails. 

The  Iterative  Method  starts  with  the  two 
curve  intersections  with  the  cell  sides.  These  are  math¬ 
ematically  connected  by  a  straight  line  segment.  To 
decide  which  additional  points  are  required,  the  line  is 
bisected  to  establish  a  reference  point  (x,y).  Then  a 
vertical  or  horizontal  line  is  constructed  through  (x,y) 
and  its  intersection  with  the  contour  is  computed.  The 
addition  of  the  new  point  to  the  original  two  means  the 
contour  can  now  be  drawn  with  two  connected  line  seg¬ 
ments.  If  the  distance  from  the  reference  point  to  the 
new  point  is  large  relative  to  the  cell's  width,  then  the 
two  new  line  segments  are  bisected  to  develop  5  points 
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along  the  curve.  New  line  segments  are  added  until  the 
distance  from  a  reference  to  its  new  point  is  "small". 
Then  the  curve  may  be  drawn  smoothly  by  connecting  the 
points  in  order  with  line  segments.  Typically,  3  to  5 
points  are  generated;  however, up  to  50  could  be  generated 
if  necessary.  The  user  influences  the  number  of  points 
generated  in  each  cell  by  defining  the  fraction  of  a  cell 
size  used  to  decide  whether  the  distance  between  the  new 
point  and  the  reference  point  is  "small". 

The  Iterative  Method  fails  when  the  algo¬ 
rithm  cannot  find  a  unique  intersection  with  the  horizon¬ 
tal  or  vertical  line  through  a  reference  point.  For  ex¬ 
ample,  if  a  curve  makes  a  Figure  "S"  across  a  cell,  then 
a  vertical  line  could  intersect  the  contour  three  times. 
The  algorithm  cannot  decide  which  to  accept.  Similar 
problems  are  encountered  when  one  contour  crosses  all 
four  sides  because  a  saddle  exists.  In  these  cases,  the 
Stepping  Method  is  used. 

The  Stepping  Method  traces  the  contour  by 
taking  very  small  steps  along  the  curve.  It  starts  at  an 
edge  point  on  the  curve  and  then  computes  elevation 
values  using  the  local  model  at  the  corners  of  sub¬ 
cells.  The  curve's  intersection  with  the  sub-cell  grid 


lines  is  computed  by  inverse  linear  interpolation.  Since 
a  curve  enters  and  exits  each  sub-cell  along  its  path,  it 
is  necessary  to  compute  local  elevation  values  for  only 
those  sub-cells.  When  the  sub-cells  are  very  small,  the 
algorithm  can  follow  the  most  complex  curve  shapes  across 
the  cell.  Sub-cells  may  be  as  small  as  1/256  the  area  of 
a  DTM  cell.  From  15  to  30  points  per  cell  are  generated 
by  this  method. 

The  Stepping  Method  is  more  time-consuming 
since  the  local  elevation  model  is  evaluated  more  times 
and  more  points  must  be  computed  across  the  cell.  It 
does  not  have  the  freedom  to  compute  only  those  points 
necessary  to  adequately  describe  the  contour.  Conse¬ 
quently,  for  efficiency  it  is  used  only  when  the  faster 
Iterative  Method  fails. 

4.4  TESTING  OF  CONTOUR-TO-GRI D  ALGORITHMS.  Testing 
of  contour-to-grid  algorithms  is  difficult  for  two 
reasons.  First,  results  must  usually  be  subjectively 
evaluated  due  to  the  lack  of  accepted  "ground  truth". 
Second,  it  is  often  difficult  to  isolate  effects  in  cases 
where  "ground  truth"  does  exist.  A  good  example  of  this 
latter  problem  could  be  created  by  the  attempt  to  evalu¬ 
ate  the  Gr id- to-Con tour  program  by  matching  the  outputs 
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of  that  program  with  the  original  contour  strings  sup¬ 
plied  by  DMA.  Any  differences  could  be  attributed  to 
either  the  Gr id- to-Contour  program  or  the  Contour- to-Gnd 
program  which  created  the  intermediate  grid.  There  is  no 
convenient  method  to  isolate  those  caused  by  the  Grid- 
to-Contour  program  alone. 

Considering  the  above  difficulties,  it  was 
decided  that  testing  would  be  performed  with  data  from 
several  small  test  areas.  These  areas  are  made  up  of  a 
small  number  of  cells  with  posted  grid  values.  For  each 
area,  the  shape  of  contours  passing  through  the  cells 
can  be  observed.  The  contours  can  be  matched  to  the  grid 
node  values  and  to  the  choice  of  intermediate  grids  over- 
layed  in  each  cell  as  is  discussed  in  Section  4.3.4.  The 
result  of  this  test  procedure  is  described  below.  How¬ 
ever,  note  that  the  contour  maps  provided  in  this  report 
represent  output  from  the  Grid- to-Contour  program.  Al¬ 
though  in  individual  cases  they  may  be  intended  primarily 
to  demonstrate  the  results  of  the  testing  performed  for 
other  parts  of  this  contract,  their  visual  satisfaction 
and  consistency  provide  added  validity  to  the  Grid-to- 
Contour  algorithms. 


The  first  test  area  shown  in  Figures  4.2 
to  4.5  are  for  small  mountainous  regions  with  many  con¬ 
tours  through  the  grid  cells.  Grid  values  at  each  node 
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:ontour  using  2x2  Overlay. 


are  posted.  In  this  case,  the  values  at  5  of  the  nodes 
correspond  to  contour  levels.  The  four  figures  corre¬ 
spond  to  4  different  intermediate  grid  choices:  2x2, 
5x5,  9x9,  and  17x17.  As  discussed  in  Appendix  B  the  2x2 
overlay  is  a  special  case  which  uses  linear  interpolation 
over  triangular  subdivisions  of  the  cell.  For  this  case 
there  are  obvious  improvements  to  be  seen  in  moving  above 
the  simpliest  2x2  overlay.  Further  improvements  are 
difficult  to  detect. 

The  second  set  of  data  shows  results  for  a 
more  difficult  area  6.  Results  are  shown  in  Figures  4.6 
to  4.9.  In  this  case,  a  small  valley  on  the  left  side 
attains  its  lowest  value  in  the  upper  left  hand  corner. 
The  higher  terrain  to  the  right  of  the  valley  protrudes 
into  the  valley  in  the  lower  left.  Again,  significant 
improvement  is  observed  in  moving  from  a  2x2  overlay  to 
higher  cases.  For  a  9x9  overlay  straight  line  segments 
in  contours  can  still  be  observed,  especially  in  the 
lower  right  corner.  For  this  particular  terrain  area  a 
17x17  overlay  is  required  to  produce  an  acceptable  out¬ 
put. 

Of  course  it  is  in  areas  with  significant 
terrain  variation  that  any  contouring  algorithm  could  be 
expected  to  perform  well.  Figures  4.10  to  4.,L3  provide  a 
more  rigorous  test  grid  lifted  out  of  test  area  4.  The 
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Figure 


4.11.  Grid- to-Contour  using  5x5  Overlay 
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Figure  4.12.  Grid- to-Contour  using  9x9  Overlay. 
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grid  is  relatively  flat  with  a  valley  deepening  and 
widening  from  left  to  right.  The  contour  at  the  height 
of  4850  presents  particular  difficulty  to  the  algorithm. 
For  this  case  the  two  coarsest  overlay  grids  supply  unac¬ 
ceptably  straight  and  sharp  shapes  to  this  contour  even 
though  the  second  level  provides  reasonable  dimensions. 
The  third  level  still  has  some  straight  contour  segments 
while  the  highest  level  finally  presents  an  acceptable 
view. 
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5.0 


SURFACE  SMOOTHING  AND  GENERALIZATION 


5.1  OVERVIEW.  The  objective  of  this  task  is  the 

development  and  testing  of  algorithms  for  smoothing 
gridded  terrain  data  and  for  generalizing  contours. 
Smoothing  may  be  defined  as  an  area  or  spatial  process 
designed  to  modify  the  significance  of  certain  types  of 
terrain  features.  Typically,  smoothing  is  performed  to 
reduce  local  surface  shaping  while  maintaining  the  re¬ 
gional  shape.  Two  major  types  of  smoothing  are  imple¬ 
mented  and  each  offers  a  number  of  options  in  applica¬ 
tion. 


Contour  generalization  is  the  process  of 
modifying  the  detail  or  placement  of  contour  lines  to 
achieve  some  display  objective.  The  operation  involves 
many  elements  and  is  ultimately  a  very  subjective  pro¬ 
cess.  This  task  is  concerned  with  only  one  element  of 
generalization,  termed  line  simplification.  This  tech¬ 
nique  is  used  when  detail  in  contour  lines  must  be  re¬ 
duced  to  accommodate  a  scale  reduction  or  to  increase  the 
clarity  of  map  features.  There  exist  a  number  ol  algo¬ 
rithms  for  the  simplification  of  lines  which  operate  by 
processing  the  line  vertices  directly.  However,  this 
task  approaches  the  simplification  of  contour  lines 
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indirectly  by  the  controlled  application  of  grid 
smoothing  prior  to  contour  line  generation. 

5.2  GRID  SMOOTHING.  Two  types  of  algorithms  have 
been  implemented  for  smoothing:  (1)  iterative  filtering 
with  a  spatial  convolution  operator  and  (2)  weighted 
least  squares  filtering.  Convolution-operator  smoothing 
is  a  recursive  process  whereby  a  small  area  operator  is 
moved  through  the  grid,  combining  values  from  both  the 
previous  pass  and  the  current  pass  to  compute  a  new  value 
at  each  grid  location.  Control  parameters  include  type 
of  operator  and  number  of  smoothing  passes.  Two  opera¬ 
tors  are  available.  These  are  a  five-point  (Laplacian) 
operator  and  a  thirteen-point  (bi-harmonic)  operator. 
The  two  operators  are  illustrated  in  Figure  5.1. 

The  Laplacian  operator  has  the  relative 
advantage  of  using  only  five  grid  values  and  thus  requir¬ 
ing  significantly  less  computation.  It  also  converges  to 
a  final  surface  very  rapidly.  Disadvantages  of  this  op¬ 
erator  are  that  the  final  grid  values  are  restricted  to 
lie  within  the  initial  elevation  limits  of  the  input  grid 
and  that  the  final  grid  has  less  than  optimal  smoothness 
properties.  By  contrast  the  bi-harmonic  operator  is 
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(a)  Laplacian  Operator 


(b)  Bi-harmonic  Operator 


Figure  5.1.  Convolution  Operators 
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derived  from  a  constraint  of  minimum  curvature  which 
guarantees  optimal  smoothness  ( sf-e  Appendix  C).  However, 
the  larger  size  and  slower  convergence  of  this  operator 
introduces  additional  calculation  expense. 

The  implementation  scheme  used  for  convo¬ 
lution  operators  is  the  technique  of  successive  over-re¬ 
laxation  (Reference  6).  The  chosen  operator  is  passed 
over  the  grid  four  times,  upper-left  to  lower  rig.:t  and 
back,  then  upper-right  to  lower-left  and  back,  to  com¬ 
plete  one  filter  cycle.  Both  newly  computed  and  previ¬ 
ously  computed  values  are  used  to  calculate  the  new  value 
at  the  current  grid  location.  The  new  value  calculated 
is  related  to  the  value  previously  calculated  by  the 
equation 


zi,j(n)  -co 


JwkZk(n-1)  +  ^Wmzm(n) 

_  K  m 


♦  ( 1-CiJ*  2  i  #  j(n-l) 


where:  n  is  the  iteration  number, 

z  is  the  gridded  value, 

Wfc  is  the  kfch  operator  weight  where  the  kfch 

location  has  recently  been  iterated, 
wm  is  the  mfch  operator  weight  where  the  m*1*1 

location  is  yet  to  oe  iterated, 

6Jis  an  acceleration  parameter. 
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Each  pass  of  the  convolution  operator 
flows  grid  information  in  the  direction  of  operator  move¬ 
ment.  To  achieve  zero  net  flow  of  the  grid  data,  a  four- 
pass  cycle  is  utilized  for  one  full  smoothing  pass. 

Weighted  least  squares  smoothing  is  pro¬ 
vided  as  a  set  of  options  that  depends  upon  local  data 
only.  Operator  shape  and  weighting  may  be  selected  from 
a  menu  of  built-in  options,  or  the  user  may  completely 
specify  the  shape  and  weighting  to  be  used.  The  only 
constraint  is  that  the  operator  must  fit  within  a  5-by-5 
grid.  The  built-in  filter  shape  options  are  shown  in 
Figure  5.2.  A  plane,  which  fits  the  grid  nodes  specified 
by  the  operator  shape  in  a  weighted  least  square  error 
sense,  is  calculated  as  described  in  Appendix  D.  The 
height  of  the  plane  at  the  location  of  the  node  to  be 
smoothed  is  used  as  the  output  value. 

5.3  CONTOUR  GENERALIZATION.  Contour  generalization 
is  used  here  in  the  restrictive  sense  of  simplifying  con¬ 
tour  lines.  It  can  be  described  as  reducing  the  amount 
of  detail  represented  by  the  contours.  The  typical  ap¬ 
proach  to  this  problem  operates  on  the  vertices  of  the 
contour  lines  directly,  producing  a  reduced  or  modified 
set  of  vertices  which  define  a  smoother  or  less  detailed 
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Figure  5.2. 
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Local  Area  Point  Specifiers  for 
Least  Square  Filtering 


curve.  There  exist  many  algorithms  of  varying  complexity 
for  defining  the  smoothing  operation  on  vertices. 

The  contour  lines  resulting  from  a  line 
smoothing  operation  may  be  thought  of  as  representing  a 
smoother  surface.  This  suggests  that  application  of  a 
grid  smoothing  technique  to  the  surface  model  would  pro¬ 
duce  similar  contours.  Simplifying  the  surface  rather 
than  the  lines  can  provide  some  advantages  such  as  pro¬ 
viding  smoothing  control  more  directly  related  to  phys¬ 
ical  attributes,  preservation  of  the  meaning  of  contour 
lines,  and  increased  efficiency  in  the  processing  of 
dense  contours. 

Any  of  the  surface  smoothing  operators 
previously  described  in  this  section  can  be  applied  for 
the  purpose  of  contour  generalization  although  additional 
control  of  the  smoothing  effect  is  necessary  for  specific 
purposes.  This  control  is  provided  in  two  forms,  fixed 
points  and  control  surfaces.  The  accommodation  of  fixed 
points  during  surface  smoothing  allows  important  features 
to  be  preserved  while  relaxing  constraints  on  the  modifi¬ 
cation  of  other  features.  Grid  nodes  which  are  identif¬ 
ied  as  fixed  points  indicate  that  the  smoothed  output 
value  for  that  point  is  to  be  identical  to  the  input 
value.  This  control  is  specified  in  the  grid  generation 
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process  described  in  Section  3.  The  points  which  define 
the  contour  lines  input  to  the  gridding  process  are  con¬ 
sidered  to  be  highly  accurate.  When  one  of  these  points 
or  the  intersection  of  a  contour  line  with  a  grid  line 
occurs  near  a  grid  node,  that  node  is  marked  as  a  fixed 
point.  The  nearness  criterion  is  usually  defined  as  half 
the  diagonal  of  a  grid  cell. 

The  second  method  of  constraining  the  sur¬ 
face  smoothing  operation  is  through  the  use  of  control 
surfaces.  These  surfaces  represent  upper  and  lower 
bounds  on  the  smoothed  surface  and  are  functions  of  the 
original  input  surface.  They  are  defined  by  their  grid 
values  as  follows: 

U(x,y)  >  S(x,y)  >  L(x,y) 

where  S(x,y)  represents  the  smoothed  surface  and  U(x,y) 
and  L(x,y)  represent  the  upper  and  lower  control  surfaces 
respectively  defined  as 

U(x,y)  =  ay  I{x,y)  +  by  , 

L(x,y)  =  aL  I(x,y)  +  bL 

where  I(x,y)  is  the  input  surface,  and  the  a's  and  b's 
are  constants  chosen  by  the  user. 
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Note  that  the  definition  of  the  control 
surfaces  does  not  require  symmetry  about  the  input  sur¬ 
face.  For  example,  peaks  may  be  restricted  while  depres¬ 
sions  are  allowed  to  fill  in.  The  two  types  of  smoothing 
constraints  are  illustrated  in  Figure  5.3. 

5.4  TESTING .  Testing  was  performed  using  the  DMA 
provided  data  of  the  Mustang  Mountains  area  in  southern 
Arizona.  As  discussed  in  Section  2  this  data  was  pro¬ 
vided  in  the  form  of  contour  strings  and  related  drain 
and  ridge  lines.  Subsets  of  this  area  with  different 
characteristics  were  chosen  as  test  areas.  These  subsets 
were  extracted  and  used  as  input  to  the  Contour- to-Grid 
program  discussed  in  Section  3.  The  output  of  this  pro¬ 
gram  is  a  65x65  grid  of  elevation  values.  The  smoothing 
program  accepts  as  input  and  produces  as  output,  grids  of 
elevation  data. 

Facets  of  the  program  tested  include  all 
major  options  of  regular  smoothing  and  constrained 
smoothing  including:  (1)  unconstrained  bi-harmonic 
smoothing,  (2)  unconstrained  Laplacian  smoothing,  (3) 
weighted  least-squares  smoothing,  (4)  constrained  smooth¬ 
ing  using  the  bi-harmonic  convolution  filter  with  control 
surfaces  or  fixed  point  constraints. 
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Figure  5.3.  Two  Methods  to  Constrain  the  Smoothing 

(a)  Control  Surfaces 

(b)  Fixed  Points 


I  1 


The  two  subset  areas  of  the  Mustang  Moun¬ 
tain  Area,  which  were  picked  for  test  of  these  algo¬ 
rithms,  illustrate  two  very  different  types  of  terrain. 
The  area  denoted  test  area  A  appears  relatively  flat  com¬ 
pared  to  the  mountainous  area  denoted  test  area  B. 

The  root  mean  square  (RMS)  curvature  sta¬ 
tistic  proved  to  be  most  useful  in  judging  the  smoothing 
performance.  For  a  DTM  with  elevation  values  Zij,  I  rows 
and  J  columns  this  statistic  is  defined  as 

J 

y  (Cij)2/i*j 
]=i 

cij  =  z i+1 ,  j  *  z  i— 1 ,  j  +  2 i  ,  j  + 1  +  z  i  ,  j  — 1  _  4Zij 

Six  contour  maps  are  shown  of  test  area  A 

in  Figures  5.4  through  5.9  and  of  test  area  B  in  Figures 
5.10  through  5.15.  They  are  arranged  in  order  of  de¬ 
creasing  curvature.  Each  smoothing  calculation  used  one 
pass. 

5.4.1  Smoothing  Constrained  by  Fixed  Points.  With 
both  data  sets  smoothing  constrained  by  fixed  points  pro¬ 
duced  only  a  small  decrease  in  RMS  curvature.  This  is 
because  the  fixed  points  represent  a  large  fraction  of 
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Figure  5.5.  Area  A  Bi-Harmonic  Smoothing  Constrained 
by  Fixed  Points 
RMS  Curvature  =  2.970 
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Figure  5.8.  Area  A  Unconstrained  Bi-Harmonic  smoothing 
RMS  Curvature  =  1.198 


71 


Area  B  Bi-Harmonic  Smoothing  Constrained 

by  Fixed  Points 

RMS  Curvature  =  13.171 
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Figure  5.13.  Area  B  Least  Squares  Smoothing 
RMS  Curvature  -  5.357 
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the  whole  grid:  27%  and  57%  for  areas  A  and  B  respec¬ 
tively. 


Comparing  both  maps  of  the  fixed  point 
constrained  smoothing  (Figures  5.5  and  5.11)  with  those 
c  the  grids  before  smoothing  (Figures  5.4  and  5.10) 
shows  only  a  slight  difference.  The  solid  contours  (at 
multiples  of  25)  remain  unchanged  since  they  are  in  re¬ 
gions  with  all  fixed  points.  This  is  because  they  follow 
the  original  contour  data  provided  by  DMA.  On  the  other 
hand  the  intermediate  contours  (dashed)  are  smoothed 
somewhat  because  these  are  away  from  the  fixed  points. 
It  appears  that  fixed  point  constraints  are  not  useful 
when  the  number  of  fixed  points  is  so  large. 

5.4.2  Least  Squares  Smoothing.  For  the  least  squares 
smoothing,  option  4  was  chosen  for  the  shape,  and  both 
weighting  functions  were  tried.  Table  5.1  gives  the  re¬ 
sulting  RMS  curvature  for  both  test  areas.  The  sharp 
weighting  gives  almost  the  entire  weight  to  the  central 
point  and  produces  almost  no  smoothing.  On  the  other 
hand,  the  smooth  weighting  produces  some  grid  averaging 
resulting  in  reduced  RMS  curvature.  The  contour  maps  of 
the  two  smooth  weighting  cases  are  shown  in  Figures  5.6 
and  5.13  They  represent  a  nice  compromise  between 
smoothness  and  detail. 
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Table  5.1 


Test  Area 

Filter 

RMS  Curvature 

A 

None 

— 

4.124 

A 

Least  Squares 

Sharp 

4.066 

A 

Least  Squares 

Smooth 

1.979 

B 

None 

— 

13.729 

B 

Least  Squares 

Sharp 

13.506 

B 

Least  Squares 

Smooth 

5.357 

5.4.3  Bi-harmonic  Smoothing  Constrained  by  Control 

Surfaces.  In  this  calculation  the  amount  of  smoothing 
produced  by  the  bi-harmonic  filter  is  limited  by  the 
control  surfaces.  For  all  such  calculations  shown  in 
this  report,  the  control  surfaces  were  set  a  constant 
distance  above  and  below  the  pre-smoothed  surface  by 
choosing  the  constants 

au  =  aL  =  i*0 
bj;  a  +b 
bL  =  -b. 

Both  test  data  sets  were  smoothed  with  many  different 
values  of  b.  Figure  5.16  shows  the  resulting  RMS  curva¬ 
ture  vs.  b.  As  b  increases,  the  RMS  curvature  decreases 
until  b=8  for  area  A  and  b=18  for  area  B,  after  which  no 
further  decrease  in  curvature  is  seen.  For  larger  values 
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RMS  Curvature 
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Figure  5.16.  RMS  Curvature  vs.  Control  Surface  Offset 
for  Con  strained  Bi -Harmonic  Smoothing 
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of  b  the  control  surfaces  no  longer  have  any  effect  and 
the  result  is  the  same  as  for  an  unconstrained  bi-harmon¬ 
ic  filter.  One  constrained  bi-harmonic  smoothing  calcu¬ 
lation  of  each  test  area  was  chosen  for  contouring.  Fig¬ 
ure  5.7  shows  the  map  of  area  A  for  b=3.0  and  Figure  5.12 
shows  the  map  of  area  B  for  b=5.0.  When  compared  to  the 
before  smoothing  maps,  these  maps  show  the  general  ter¬ 
rain  features  much  more  clearly. 

5.4.4  Bi-harmonic  Smoothing.  When  not  constrained, 
the  bi-harmonic  filter  produces  a  large  decrease  in  RMS 
curvature,  yet  the  important  terrain  features  are  pre¬ 
served.  Figures  5.8  and  5.14  show  the  appropriate  con¬ 
tour  maps.  This  filter  acts  quite  fast  for  three  rea¬ 
sons:  (1)  It  is  a  recursive  filter  acting  on  both  old  and 
new  grid  values,  (2)  Each  pass  through  the  filter  in  fact 
represents  four  sub-passes  (necessary  to  remove  any  phase 
shift),  and  (3)  The  acceleration  parameter  (1.3)  was 
chosen  for  high  speed. 

5.4.5  Laplacian  Smoothing.  Use  of  the  Laplacian 
smoothing  filter.  Figures  5.9  and  5.15  results  .n  a  much 
more  drastic  smoothing  effect  for  one  pass.  The  terrain 
details  are  thoroughly  washed  out,  and  only  the  general 
trends  remain.  Since  the  Laplacian  operator  is  smaller 
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than  the  bi-harmonic  operator,  it  requires  less  computa¬ 
tion,  and  it  may  be  preferred  for  its  great  amount  of 
smoothing . 

5.5  CUTOFF  OF  THE  BI-HARMONIC  OPERATOR.  The  bi¬ 

harmonic  smoothing  filter  is  often  iterated  through 
several  passes.  This  may  become  very  time  consuming  for 
large  grid  size.  Four  statistics  were  examined  to  find  a 
criterion  to  limit  the  number  of  passes.  The  statistics 
are: 

1)  RMS  Curvature 

2)  Change  in  RMS  Curvature 

3)  Standard  Deviation 

4)  ^Z  max  (i.e.  largest  change  in  grid 

value ) 

5)  f\z  max/Range  (i.e./^Z  divided  by  grid 

maximum-grid  minimum) 

Table  5.2  and  5.3  tabulate  these  five  quantities  for 
areas  A  and  B  respectively.  These  are  the  results  after 
each  fourth  subpass,  since  it  is  always  required  to  fil¬ 
ter  in  multiples  of  four  subpasses  to  avoid  introducing 
distortions.  See  Appendix  C. 

For  all  four  constrained  cases,  most  of 
the  change  occurred  during  the  first  pass.  The  standard 


Table  5.2. 

Statistics  for  Bi-Harmonic  Smoothing  of  Area  A 


Constraint 

Pass 

RMS 

Curvature 

Change  in 

RMS 

Curvature 

Standard 

Deviation 

Az  max 

A  Z  max 

Range 

0 

4. 124 

'  ' 

25.00 

' 

Fixed 

1 

2.970 

28.0% 

25.04 

2.678 

.01847 

Point 

2 

2.924 

1.5% 

25.07 

.787 

.00559 

3 

2.914 

.3% 

25.08 

.392 

.00280 

1 

4 

2.909 

.2% 

25.10 

.278 

.00199 

1 

5 

2.907 

.1% 

25.11 

.226 

.00162 

1 

L 

6 

2.905 

.1% 

25. 12 

.183 

.00131 

1 

! 

0 

4.124 

1 

25.00 

|  Control 

1 

1.971 

52.2% 

24.93 

1.215 

.00838 

|  Surface 

‘2 

1.921 

2.5% 

24.93 

.383 

.00264 

|  b-5.0 

3 

1.914 

.4% 

24.92 

.271 

.00187 

1 

4 

1.911 

.2% 

24.92 

.144 

.00099 

1 

5 

1 .909 

.1% 

24.91 

.  106 

.00073 

1 

l 

6 

1.908 

.  1% 

24.90 

.091 

.00063 

1 

1 

' 

0 

4.124 

25.00 

|  None 

1 

1 .  198 

71  .0% 

24.87 

2.738 

.01888 

1 

2 

.873 

27.1% 

24.82 

1.095 

.00768 

1 

3 

.761 

12.8% 

24.79 

.654 

.00463 

1 

4 

.693 

8.9% 

24.76 

.498 

.00355 

1 

5 

.645 

6.9% 

24.73 

.398 

.00285 

1 

i _ 

6 

1 _ 

.609 

_ 

5.6% 

24.71 

.325 

.00233 

Table  5.3 


Statistics  for  Bi-Harmonic  Smoothing  of  Area  B 


Constraint 

Pass 

RMS 

Curvature 

Change  in 

RMS 

Curvature 

Standard 

Deviation 

A  Z  max 

AZ  max 

Range 

'  ' 

0 

13.729 

161.04 

Fixed 

1 

13. 171 

4.  i« 

161.15 

7.414 

.00925 

Point 

2 

13. 151 

0.2% 

161.14 

1.459 

.00182 

3 

13.147 

0.0% 

161.13 

.563 

.00070 

4 

13.146 

0.0% 

161.12 

.334 

.00042 

5 

13.145 

0.0% 

161.11 

.248 

.00031 

6 

13. 145 

0.0% 

. 

161.11 

.188 

.00023 

0 

13.729 

161.04 

Control 

1 

7.329 

46.6% 

161.02 

2.309 

.00288 

Surface 

2 

7.286 

0.6% 

161.03 

.577 

.00072 

b-5.0 

3 

7.280 

0.1% 

161.04 

.315 

.00039 

4 

7.278 

0.0% 

161.05 

.253 

.00032 

5 

7.277 

0.0% 

161.06 

.226 

.00028 

6 

7.276 

. 

0.0% 

161.07 

.  195 

.00024 

0 

13.729 

161.04 

None 

1 

3.526 

74.3% 

160.98 

9.641 

.01203 

2 

2.371 

32.8% 

160.93 

2.354 

.00294  | 

3 

2.110 

11.0% 

160.89 

1.422 

.00177  j 

4 

1.968 

6.7% 

160.86 

1.052 

.00131  j 

5 

1.864 

5.3% 

160.84 

.848 

.00106  | 

6 

1.782 

4.4% 

160.81 

.715 

.00089  | 
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deviation  measurement  changes  only  slightly  and  does  not 
decrease  uniformly,  making  it  useless  as  ».  cutoff  crite¬ 
rion.  The  change  in  RMS  curvature  and  the  Az  parameters 
indicate  that  two  or,  at  most,  three  passes  may  be  use¬ 
ful.  The  RMS  curvature  seems  to  be  the  most  relevant 
measure  since  this  is  the  value  which  the  filter  seeks  to 
minimize.  Selecting  a  cutoff  value  between  .5  percent 
and  2  percent  would  suffice  for  all  constrained  cases 
tabulated.  Using  the  much  less  costly  calculation  of  AZ 
Max/Range,  a  cutoff  value  of  approximately  .003  would 
give  similar  results. 

The  unconstrained  cases  are  interesting 
since  they  may  indicate  how  the  filter  is  behaving  in  the 
free  areas  of  the  constrained  grid.  Again  the  standard 
deviation  is  useless.  The  change  in  RMS  curvature  would 
provide  the  best  control  and  appears  consistent  for  the 
two  cases.  Although  the  AZ  Max/Range  is  decreasing 
steadily,  a  fixed  cutoff  value  would  give  a  very  differ¬ 
ent  performance  for  the  two  unconstrained  cases.  Over¬ 
all,  the  change  in  RMS  curvature  seems  to  provide  the 
best  indicator  of  filter  effects.  However,  one  caution 
is  necessary  in  the  use  of  this  measure.  Since  the 
calculation  is  an  average  of  curvature  over  the  entire 
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grid,  significant  changes  in  a  small  area  might  be 
masked.  At  the  other  extreme  the  AZ  Max/Range  criterion 
measures  the  change  at  only  one  grid  node.  For  the 
constrained  cases,  the  behavior  of  these  two  indicators 
seems  to  correspond  well  enough  to  suggest  that  either 
would  provide  adequate  results.  Since  the  AZ  Max/Range 
requires  much  less  computation,  it  is  recommended  for  use 
as  the  filter  cutoff  criterion. 
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6. 


RESAMPLING  STUDY 


6.1  OVERVIEW.  This  part  of  this  study  compares  the 

performance  of  two  interpolation  algorithms  in  the  resam¬ 
pling  of  digitized  terrain  models.  The  algorithms  inves¬ 
tigated  are  the  local  bicubic  surface  fitting  method  de¬ 
veloped  by  Akirna  (Reference  3  and  4)  and  the  ZYCOR  para¬ 
bolic  adaption  to  the  Jancaitis  method  (Reference  2). 
Detailed  descriptions  of  these  algorithms  appear  in  Ap¬ 
pendices  F  and  G.  Comparison  between  the  two  methods  was 
based  primarily  on  the  statistical  performance  of  the 
algorithms  in  resampling  terrain  data  provided  by  DMA 
representing  the  Mustang  Mountain  area  in  southern  Ari¬ 
zona  and  visual  inspection  of  results.  Also  compared  are 
the  spatial  frequency  characteristics  of  the  algorithms, 
their  response  to  impulses,  and  the  amount  of  time  re¬ 
quired  by  each  to  perform  interpolations. 

Four  basic  conclusions  are  available  from 

this  study: 

1)  The  two  algorithms  provide  very  simi¬ 
lar  performance  in  resampling  the  test 
terrain  provided  by  DMA.  While  local 
differences  where  observable  in  the 
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resampling  of  individual  DTMs,  these 
did  not  point  to  overall  patterns 
which  would  lead  to  choosing  either 
algorithm  as  generally  superior  to  the 
other . 

2)  The  Jancaitis  algorithm  runs  signifi¬ 
cantly  faster  than  the  Akima  algorithm 
as  implemented.  Based  on  test  cases 
with  input  grids  generated  by  random 
number  routines  and  running  on  ZYCOR's 
VAX  11/750,  the  Jancaitis  algorithm 
out-performed  the  Akima  routine  at  a 
rate  of  2.5  to  1. 

3)  The  Akima  algorithm  can  perform  a 
transition  into  a  terrain  with  zero 
curvature  without  the  overshoot  pro¬ 
blem  produced  in  the  use  of  the 
Jancaitis  algorithm.  A  study  of  the 
response  of  both  algorithms  to  a  sin¬ 
gle  impulse  provides  a  quantification 
of  this  difference. 

4)  A  comparison  of  the  power  spectra  of 
the  data  grids  and  that  of  grids  pro¬ 
duced  by  resampling  at  the  same  den¬ 
sity  using  the  two  interpolation 
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algorithms  showed  a  small  advantage 
for  the  Jancaitis  method  as  measured 


by  the  attenuation  of  high  frequencies. 
All  of  the  above  points  are  elaborated  on  the  following 
sections. 

6.2  RESAMPLING  OF  DMA  DATA.  As  discussed  in  Section 
2,  DMA/ETL  provided  ZYCOR  with  digitized  contour  data  for 
7.5  minute  x  7.5  minute  area  in  southern  Arizona  around 
Mustang  Mountain  near  Fort  Huachuca.  This  data  was  parti¬ 
tioned  by  ZYCOR  into  a  number  of  smaller  test  areas  repre¬ 
senting  various  terrain  types  and  the  partitioned  areas 
were  then  gridded  to  provide  test  DTMs .  Table  6.1  gives 
the  grid  sizes  and  gridding  increments  in  both  coordinates 
for  the  four  test  areas  considered  in  this  section.  The 
numbering  scheme  which  runs  1,4,5,  and  6  reflects  the  fact 
that  certain  areas  originally  chosen  for  testing  were 
eventually  disgarded  due  to  difficulties  in  running  CTOG 
on  the  data  or  problems  with  the  data  itself.  These  prob¬ 
lems  have  been  identified  as  mistakes  made  in  the  use  of  a 
program  to  partition  the  data  sets. 

In  order  to  provide  statistical  measures 
of  accuracy  it  was  desired  to  resample  the  grids  describ¬ 
ed  above  at  the  actual  grid  node  locations.  To  do  this 


90 


Table  6.1. 

Test  Grids  Used  in  Resampling  Study 


Test  Area  1 

Grid  Size: 

48  x  48 

Gridding  Increments: 

78m  x  92m 

Test  Area  4 

Grid  Size: 

48  x  48 

Gridding  Increments: 

52m  x  61m 

Test  Area 

5 

Grid  Size: 

64  x  64 

Gridding  Increments: 

58m  x  68m 

Test  Area  6 

Grid  Size: 

48  x  48 

Gridding  Increments: 

58m  x  68m 

each  grid  node  was  assumed  to  be  the  center  point  of  a 
cell  in  a  thinned  grid  created  by  the  ignoring  of  every 
other  row  and  column  of  the  original  including  the  row 
and  column  which  define  the  grid  node.  Figure  6.1  shows 
this  sampling  scheme.  In  order  to  avoid  possible  non¬ 
general  effects  of  using  these  algorithms  to  interpolate 
near  edges  of  a  grid,  resampling  was  restricted  to  take 
place  in  the  center  portion  of  a  grid  ignoring  the  five 
rows  and  columns  along  each  border. 

Considering  the  fact  that  the  center  of  a 
cell  may  not  be  the  point  at  which  the  most  significant 
differences  between  the  two  resampling  algorithms  may  be 
observable  (see  below,  Section  6.4),  a  second  method  was 
developed  which  involved  retaining  as  data  only  every 
third  row  and  column  of  the  input  grid  as  data.  This 
resulted  in  12  possible  interpolation  locations  relative 
to  a  grid  cell.  However,  results  obtained  with  this 
method  in  no  way  modified  or  added  to  results  obtained  in 
using  the  method  described  above,  so  this  approach  will 
not  be  pursued  further  in  this  report. 

During  the  running  of  the  resampling  pro¬ 
gram  the  difference  between  the  interpolated  value  at 
each  grid  node  and  the  actual  value  was  computed.  The 
errors  were  used  to  compute  three  statistical  measures  of 
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performance  by  which  the  two  interpolation  algorithms 


could  be  compared:  the  square  root  of  the  average  of  the 
squared  errors  (RMS),  the  average  of  the  absolute  errors, 
and  the  error  with  the  largest  absolute  magnitude.  The 
interpolated  grids  were  output  in  a  formatted  form  so 
that  they  could  be  contoured  and  plotted  for  comparison 
against  the  original  data.  Furthermore,  the  difference 
grids  between  the  interpolated  and  original  grids  were 
output  for  both  the  Akiina  and  Jancaitis  algorithms.  Also 
output  were  the  norm  of  the  gradient  and  the  absolute 
curvature  for  each  node  of  the  input  grid. 

For  grid  node  position  i,j  the  absolute 
value  of  the  curvature  is  deterred  by 

|  zi+-l']  +  2  i-l'j  +  zi'3  +  l  +  4  Zil  | 

and  the  norm  of  the  gradient  by 


1  i  +  1 > ]  z  l-l  '  j 

2  + 

zi']  +  l  "  z  i '  j-l 

_  2  — 

—  2  — 

The  process  of  throwing  out  rows  and  col¬ 
umns  during  the  resampling  forced  our  attention  to  the 
spatial  frequency  characteristics  of  the  data.  If  a  sig¬ 
nificant  portion  of  the  inherent  frequencies  of  the  input 
grid  were  too  high,  then  no  possible  interpolation  jheme 
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which  involved  ignoring  some  of  the  input  data  was  going 
to  provide  reasonable  results.  We  therefore  determined 
to  eliminate  these  high  frequency  terms  from  the  data. 
Since  the  test  scheme  involved  throwing  out  every  other 
row  and  column,  it  was  judged  that  the  high  frequency 
terms  needing  elimination  were  those  above  fs/2  where  fs 
was  the  foldover  frequency  implied  by  the  sampling  rate 
of  the  original  grid.  A  Digital  Frequency  Transform 
(  DFT )  was  provided  to  test  and  smooth  the  input  data. 
The  transform  used  was  based  on  a  well-known  algorithm 
for  carrying  out  Fast  Fourier  Trans  forms ( FFT )  (Reference 
5.)  For  each  input  grid  a  DFT  was  performed  on  each  row. 
High  frequency  terms  which  would  be  aliased  by  the  sampl¬ 
ing  scheme  were  set  to  zero.  The  spectrum  for  each  row 
was  then  inverse  transformed  back  to  the  spatial  domain. 
This  procedure  was  then  repeated  on  the  columns.  This 
process  produced  a  smoother  grid  with  the  most  noticeaole 
difference  being  radical  changes  around  the  edges  of  the 
grid.  These  differences  were  attributable  to  the  cycli¬ 
cal  nature  of  performing  frequency  analysis  with  sampled 
data  which  leads  some  of  the  high  frequency  components  of 
the  spectrum  to  result  from  differences  between  the  end 
of  each  row  and  column  of  the  grids.  By  throwing  out 
high  frequencies  these  ends  were  forceably  adjusted  to 
fit  each  other. 
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A  summary  of  the  results  of  these  test 
efforts  is  contained  in  Table  6.2.  Contour  plots  corre¬ 
sponding  to  these  tables  are  contained  in  Figures  6.2 
through  6.33.  For  each  test  area,  eight  grids  are  shown. 
These  contain: 

1)  the  original  test  area, 

2)  the  DFT  smoothed  test  area, 

3)  the  smoothed  test  area  resampled 

using  the  Akima  algorithm, 

4)  the  smoothed  test  area  resampled  us¬ 
ing  the  Jancaitis  algorithm, 

5)  an  absolute  error  grid  for  the  Akima 
algorithm, 

6)  an  absolute  error  grid  for  the 

Jancaitis  algorithm, 

7)  a  contoured  grid  of  the  norm  of  the 

smoothed  grid  gradient,  and 

8)  a  contour  grid  of  the  absolute  curva¬ 
ture  of  the  smoothed  grid. 

The  number  contained  in  Table  6.2  shows 

that  little  quantitative  difference  was  produced  from 
this  testing  procedure.  In  all  but  one  case  the  maximum 
difference  in  the  error  statistics  produced  by  the  two 
methods  was  less  than  10  percent.  In  most  cases  they 
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Table  6.2 

Grid  Resampling  Statistics 


Test  Area  1 

ERRORS 

AKIMA 

JANCAITI S 

RMS 

+  3.96 

3.67 

AVR 

2.44 

2.30 

PEAK 

-35.00 

-33.9 

Test  Area  4 

ERRORS 

AKIMA 

JANCAITI S 

RMS 

4.89 

4.67 

AVR 

3.40 

3.23 

PEAK 

-31.0 

-33.4 

Test  Area  5 

ERRORS 

AKIMA 

JANCAITI S 

RMS 

10.3 

10.5 

AVR 

7.1 

7.16 

PEAK 

-52.1 

-78.4 

Test  Area  6 

ERRORS 

AKIMA 

JANCAITIS 

RMS 

17.6 

17.8 

AVR 

13.6 

13.5 

PEAK 

-79. 

-81.0 
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6.2.  Test  Area  1 
5  Foot  Contours 


Figure  6.4.  Akima  Interpolation  of  Test  Area  1 
25  Foot  Contours 
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Figure  6.7.  Jancaitis  Error  Grid  for  Test  Area  1 

5  Foot  Contours 


Figure  6.8.  Absolute  Curvature  Test  Area  1 
20  Foot  Contour* 
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6.10.  Test  Area  4 
Foot  Contours 
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.  DFT  Smoothed  Test  Area  4 
25  Foot  Contours 
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.12.  Akima  Interpolation  of  Test  Area  4 
25  Foot  Contours 
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Figure  6.15.  Janc'.itis  Error  Grid  Test  Area  4 
5  Foot  Contours 
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Figure  6,16.  Absolute  Local  Curvature  Test  Area  4 

20  Foot  Contours 
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Norm  of  Surface  Gradient  Test  Area  4 
10  Foot  Contours 


Figure  6.20.  Akima  Interpolation  Test  Area  5 
50  Foot  Contours 
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Figure  6.22.  Akima  Error  Grid  Test  Area  5 
10  Foot  Contours 
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Figure  6.28.  Akima  interpolation  Test  Area  6 
50  Foot  Contours 
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Figure  6.30.  Akima  Error  Grid  for  Test  Area  6 
10  Foot  Contours 
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were  less  than  5  percent.  The  one  example  of  an  error 
greater  than  10  percent  is  a  measurement  of  peak  error 
for  which  the  largest  variations  are  to  be  expected. 
Visual  differences  between  the  results  for  the  two 
algorithms  were  also  very  small.  Without  careful 
observation  of  the  interpolated  grids  or  the  generated 
error  grids,  it  is  not  easy  to  distinguish  differences 
between  the  two  algorithms. 

At  first  it  may  seem  surprising  that  two 
algorithms  which  are  so  different  in  derivation  would  not 
produce  more  noticeably  different  results.  However, 
although  there  are  great  differences  in  implementation  of 
the  two  algorithms,  there  are  certain  basic  conditions 
which  make  the  similarity  of  results  reasonable.  These 
are : 


1) 

Both  algorithms 

are  local 

operators. 

2  ) 

The  interpolated 

surface  generated  by 

both  algorithms 

must  match 

the  origi- 

nal  grid  at  corners. 

3) 

The  transitions 

between 

cells  for 

both  algorithms 

must  be 

smooth  as 

measured  by  first  derivatives. 

4)  The  larger  input  area  used  in  carry¬ 
ing  out  the  Akima  algorithm  is 
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balanced  against  the  higher  degree 
polynomial  used  in  the  Jancaitis 
algorithm. 

6.3  TIMING  CONSIDERATIONS.  The  relatively  small 
performance  differences  observed  between  these  two  algo¬ 
rithms  lead  naturally  to  attempting  to  differentiate 
between  them  on  the  basis  of  their  speed  of  execution. 
To  this  end,  two  small  test  programs  were  written  which 
generated  test  grids  by  use  of  a  random  number  generator 
and  then  called  the  individual  interpolation  algorithms 
repeatedly.  The  Jancaitis  algorithm  represents  code 
entirely  generated  by  2YC0H  and  optimized  to  provide 
rapid  computations  through  re-use  of  interpolation 
weights.  The  Akima  algorithm  used  in  these  tests  was 
based  on  published  code.  In  order  to  make  the  comparison 
between  the  two  reasonable,  the  Akima  code  was  modified 
to  make  it  a  function  instead  of  a  subroutine  and  it  was 
made  applicable  only  to  unit  grids.  This  last  modifica¬ 
tion  removed  multiple  divisions  existing  in  the  original 
code  which  are  not  required  in  grids  with  identical  di¬ 
mensions  in  both  directions. 

In  order  to  prevent  the  exact  form  of  the 
grid  generated  by  the  random  number  generator  from  af¬ 
fecting  the  results,  the  tests  were  run  in  sections  of 
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100  interpolations  with  a  new  grid  being  created  after 
each  100  interpolations  and  the  original  seed  number  for 
the  random  number  generator  being  created  by  a  built  in 
VAX  function  which  returns  clock  time  to  a  calling  pro¬ 
gram. 

A  run  of  5000  total  interpolations  re¬ 
quired  21.6  CPU  seconds  on  ZYCOR's  VAX  11/750  with  the 
Akima  algorithm  and  8.18  CPU  seconds  with  the  Jancaitis 
algorithm.  The  advantage  of  the  Jancaitis  algorithm  is 
to  be  expected  from  a  comparison  of  the  number  of  float¬ 
ing  point  operations  required  by  the  two  algorithms.  The 
Jancaitis  subroutine  uses  56  floating  multiplications  and 
divisions  and  64  floating  point  additions  while  the  Akima 
algorithms  uses  78  multiplications  and  139  divisions.  A 
hardware  floating  point  accelerator,  not  available  for 
the  VAX  11/750  at  the  time  of  these  tests,  might  reduce 
the  ratio  of  run  times  by  a  small  amount. 

It  should  be  noted  that  much  of  the  extra 
time  required  by  the  Akima  algorithm  is  used  in  computing 
the  required  partial  derivatives.  In  applications  where 
a  DTM  is  being  sampled  at  a  rate  at  least  as  fine  as  its 
own  grid  and  in  which  the  interpolation  points  can  be 
ordered  much  time  could  be  saved  by  saving  and  reusing 
those  partial  derivatives.  This  is  discussed  further  in 
Appendix  B. 
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6.4  TRANSITION  INTO  FLAT  TERRAIN.  The  Akima  algo¬ 
rithm  can  move  from  rough  terrain  into  flat  terrain  with¬ 
out  overshoot.  This  is  an  advantage  not  displayed  by  the 
Jancaitis  algorithm.  Figure  6.34  and  the  following  dis¬ 
cussion  quantify  this  difference.  The  figure  shows  the 
result  of  applying  the  two  algorithms  to  a  profile  of  a 
grid  consisting  of  0's  except  for  a  single  isolated  grid 
value  of  height  D  (i.e.  a  D  unit  impulse).  The  curved 
surfaces  represent  the  values  produced  using  the  two 
algorithms  to  interpolate  to  all  intermediate  points 
between  the  grid  nodes  along  the  profile.  The  polynomial 
function  defining  the  surface  between  each  pair  of  nodes 
is  also  indicated.  ' 

The  top  part  of  the  figure  shows  that  the 
Akima  algorithm  fits  a  smooth  curve  out  to  the  first  0 
grid  node  and  then  is  exactly  zero  from  then  on.  The 
first  derivative  of  the  interpolation  polynomial  is  0  at 
the  top  of  the  impulse  and  at  the  grid  nodes  on  either 
side  of  it.  This  rapid  pickup  of  the  flat  terrain  is  due 
to  the  intelligent  way  the  weights  are  applied  to  the 
pair  of  different  equations  used  to  estimate  the  slope  at 
individual  nodes  (see  Appendix  B).  The  fact  that  the 


surface  around  node  2  is  flat  with  0  curvature,  implies 
that  no  weight  is  applied  to  the  difference  between  grid 


A  k  l  m  a  Interpolation  Impulse  kespcnse 


,  d*(-3*ixl5  +  7.5*|xl4-4.5*lxi^-lxl2  +  l) 


Jancai  t i s  Interpolation  Impulse  Res  po  n  s  e 


Figure  fa. 14 


Interpolation  Overshoot 
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values  at  locations  0  and  1  and  that  full  weight  is  ap¬ 
plied  to  the  0  difference  between  grid  nodes  1  and  2. 

In  the  lower  part  of  the  figure  we  see 
that  the  Jancaitis  algorithm  fits  a  smooth  curve  out  to 
the  first  0  grid  node  and  then  dips  below  the  0  level  to 

achieve  a  minimum  value  of  -1  D  at  1  +  1/  y/To. 

12 

The  slope  of  the  interpolation  surface  at  grid  location  1 
is  -d/2.  The  Jancaitis  algorithm  produces  overshoot 
throughout  this  interval. 

This  difference  between  the  two  algorithms 
will  repeat  in  a  slightly  modified  form  for  any  transi¬ 
tion  from  a  rough  part  of  a  grid  into  one  with  zero  cur¬ 
vature. 

6.5  FREQUENCY  POWER  SPECTRUM  STUDY.  In  order  to 
focus  on  frequency  questions  without  regard  to  aliasing, 
etc.,  it  was  determined  to  do  more  tests  by  resampling 
the  grids  at  the  center  of  each  cell  using  both  algo¬ 
rithms.  The  result  was  a  resampled  version  of  the  input 
grid  for  which  any  lost  information  would  be  attributable 
only  to  the  interpolation  schemes  themselves. 

Power  spectra  for  the  original  grids  and 
the  grids  produced  by  the  interpolation  methods  were  com¬ 
puted  and  compared.  To  compute  the  spectra,  the  follow¬ 
ing  steps  were  followed: 


1)  Each  individual  row  and  column  was 
normalized  so  that  its  mean  was  0  and 
its  standard  deviation  was  1.  This 
step  was  intended  to  prevent  a  few 
rows  or  columns  from  dominating  the 
final  results. 

2)  Each  row  and  column  was  tapered  using 
a  Hamming  window  to  reduce  sidelobes 
created  by  taking  a  finite  sample  of 
terrain  data  and  by  the  cyclical  na¬ 
ture  of  a  DFT. 

3)  A  Digital  Frequency  Transform  of  each 
row  and  column  was  computed  and  the 
in-phase  and  quadrature  terms  squared 
and  added  to  produce  a  power  spectrum 
for  the  row  or  column. 

4)  The  power  spectra  for  all  rows  and 
columns  were  added  to  provide  an  com¬ 
bined  power  spectrum  for  the  entire 
grid. 

Figures  6.35  through  6.42  show  the  resam¬ 
pled  grids  produced  for  this  part  of  the  study.  Tables 
6.3  to  6.6  show  the  power  spectra  of  the  original  sur¬ 
face,  the  Akima  resampled  surface,  the  Jancaitis  resam¬ 
pled  surface,  and  the  ratio's  of  the  Akima  and  Jancaitis 
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Figure  6.35 


Akima  Resampling  of  Test  Area  1  for 
Frequency  Analysis 
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re  6.38.  Jancaitis  Resampling  of  Test  Area  4  for 
Frequency  Analysis 
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Figure  6.40  Jancaitis  Resampling  of  Test  Area  5  for 

Frequency  Analysis 
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.42.  Jancaitis  Resampling  of  Test  Area  6 
Frequency  Analysis 
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Table  6  .  3 


Power  Spectrum  Test  Results  for  Test  Area  4 
Length  of  Power  Spectrum  *  20 


F'OUF K  SFECTRUH  RESULTS 


N 

RAM  DATA  GRID 

AKIMA  GRID 

JANCAI T IS  GRID 

A,  R 

1 

1510.400 

1520.400 

1512.700 

1.007 

n 

5719.900 

5873.200 

5847.800 

1.027 

3 

2217.400 

2347.000 

2328.300 

1.058 

4 

1007.300 

1015.500 

1029.500 

1.008 

5 

730.570 

755.010 

772 .950 

1.033 

6 

892.560 

961 .800 

967.510 

1.078 

7 

568.310 

580.110 

592.220 

1 .021 

8 

389.680 

341 .440 

365.180 

0.876 

9 

367.840 

286.350 

312.650 

0.778 

10 

210.180 

159.700 

175.130 

0.760 

11 

125.020 

81.810 

88.925 

0.654 

12 

101.750 

54.157 

58.452 

0.532 

13 

84.692 

37.527 

40.376 

0.443 

14 

75.265 

21.306 

24.171 

0.283 

15 

56.157 

11  .6R5 

13.194 

0.208 

16 

45.290 

8.875 

9.203 

0.  196 

17 

50.534 

6.193 

6.05? 

0.123 

18 

42 .  '*25 

3.407 

2.740 

0.081 

19 

30.’  667 

2.561 

1.441 

0.0B4 

20 

27.947 

2.303 

0.918 

0.079 

J/R 

1.00? 
1.022 
1.020 
1.022 
1.0GB 
1 .084 
1.04? 
0.93? 
0.8G0 
0.833 
0.711 
0.G74 
0.477 
0.321 
0.23G 
0.203 
0.120 
0.065 
0.047 
0.033 
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Table  6 . 4 


Power  Spectrum  Test  Results  for  Test  Area  1 
Length  of  Power  Spectrum  =  28 


POWER  SPECTRUM  RESULTS 


N 

RAW  DATA  GRID 

AKIMA  GRID 

JANCAITIS  GRID 

A/R 

J/R 

1 

1652.100 

1443.500 

1441 .900 

0.874 

0.873 

17232.000 

15807.000 

15758.000 

0.917 

0.914 

3 

3709.600 

3949.400 

3942.000 

1.065 

1.063 

4 

1842.600 

1  *758 . 000 

1966.100 

1.063 

1.067 

C 

■J 

809.290 

820.780 

825.700 

1.014 

1 .020 

6 

350.730 

337 . 780 

341 .060 

0.962 

0.972 

7 

292.650 

276.810 

282.790 

0.946 

0.966 

8 

149.100 

123.010 

128.310 

0.825 

0.861 

9 

187.520 

160.980 

168.990 

0.858 

0.901 

10 

145.840 

131,450 

137.770 

0.901 

0.945 

1 1 

99.788 

75.131 

79.783 

0.753 

0.800 

12 

131.160 

117.150 

125.550 

0.893 

0 . 957 

13 

137.440 

110.360 

118.680 

0.803 

0.864 

14 

71 . 267 

49.094 

54.555 

0.689 

0 . 766 

15 

53.821 

30.272 

74,120 

0.562 

0.634 

16 

58 . 667 

32.658 

36  *  655 

0.557 

0.625 

1  7 

31.472 

14.534 

16.420 

0.462 

0.528 

18 

21 .299 

8.469 

9.301 

0.398 

0.437 

19 

22.682 

7.463 

7.624 

0.329 

0.336 

70 

34.837 

10.274 

11.979 

0.295 

0.344 

71 

36 . 825 

6.307 

7.840 

0.235 

0.292 

i  -) 

17.669 

3.167 

3.689 

0.179 

0.209 

23 

19.370 

4  .  1  42 

4.719 

0.214 

0.244 

2  4 

14.387 

1.988 

2 . 205 

0.138 

0.153 

25 

14.755 

1.826 

1 .640 

0.124 

0.111 

76 

15.651 

1.141 

1.244 

0.073 

0.079 

27 

14.501 

1.343 

1  .123 

0.093 

0.077 

28 

15.033 

1.264 

0.872 

0.084 

0.058 
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Table  6.5 


Power  Spectrum  Test  Results  for  Test  Area 
Length  of  Power  Spectrum  =  28 


POWER  SPECTRUM  RESULTS 


W  DATA  GRID 

AKIMA  GRID 

JANCAIT1S  GRID 

A/R 

7080.500 

7050.100 

7009.900 

0.996 

20319.000 

20411.000 

20342.000 

1 .006 

4522.800 

4870.200 

4765.800 

1.07/ 

2261.700 

2577.100 

2503.200 

1.13V 

1136.800 

1216.400 

1195,400 

1  .070 

527.420 

532.280 

522.700 

1 .009 

437.930 

462.600 

458.150 

1 . 056 

575.700 

600.240 

600.100 

1  .043 

554.290 

615.850 

613.100 

1.111 

350.980 

363.050 

361.020 

1.034 

189.790 

178.820 

182.280 

0.94? 

228.470 

197.560 

212.110 

0.865 

318.950 

264.740 

284.830 

0.830 

275.530 

222.240 

236.660 

0.807 

177.490 

135.370 

145.130 

0.763 

124.400 

86.062 

92.880 

0.692 

99.102 

47.685 

52.830 

0.481 

89.793 

34.180 

38.907 

0.384 

81.520 

25.381 

28.207 

0.311 

73.359 

23.128 

26.233 

0.315 

71.904 

19 . 441 

21.201 

0.270 

62.127 

10.507 

12.148 

0.169 

40.334 

6.830 

6.937 

0.169 

35.460 

5.324 

4.905 

0.150 

45.162 

4.307 

4.382 

0.095 

37.864 

3.591 

2.559 

0.095 

35.006 

3.233 

1.753 

0.092 

39.535 

3.407 

1,670 

0.086 

J/R 

0.990 
1.001 
1.054 
1.107 
1  .052 
0.991 
1.044 
1.042 
1.106 
1.029 
0.960 
0.930 
0.893 
0.859 
0.818 
0.747 
0.533 
0.433 
0.346 
0.350 
0.295 
0.196 
0.17? 
0.138 
0.097 
0.068 
0.050 
0.042 
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Table  6.6 


Power  Spectrum  Test  Results  for 
Lenyth  of  Power  Spectrum 


Test  Area  6 
=  20 


F'UWfft  SPEC  fRUtt  RESULTS 


N 


1 

o 

3 

4 

5 

6 
1 
8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


PAW  DATA  GRID  AMhA  GRID  JANCMTIS  GRID 


1341.900 
5788.300 

2080.900 
606.760 
281.980 
117.930 

52.993 

61.427 

34.921 

70.382 

13.938 

13.131 

7.206 

6.066 

5.928 

4.292 

5.117 

3.895 

3.568 

3.925 


14  76 
5828 
2078 
576 
265 
107 
48 
55 . 
27, 
16, 
8. 
7. 
3. 
3. 
2. 
-I 

1  .’ 
1  . 
1  . 


.  900 
.600 
.  400 
.  350 
.850 
.410 
.350 
.639 
.130 
.  710 
.854 
.649 
.915 
.22  7 
’433 
,126 
637 
,555 
217 
126 


1479.500 
583 1 . 700 
2081 . V00 
577.330 
269,720 
109.870 
48.795 
*‘6.665 
.7.585 
17.411 
8.991 
7.864 
4.06? 
3.094 
2.308 
1,958 
1.619 
1.319 
1.135 
1.072 


ft /ft 

1.101 

1  .007 

0.999 

0.950 

0,943 

0.911 

0.91? 

0.864 

0.777 

0.870 

0.635 

0.582 

0.543 

0.532 

0.409 

0.495 

0,319 

0.399 

0.341 

0.287 


J/R 


1  .103 
1.007 
1 .000 
0.951 
0.955 
0.932 

0.9?i 

0.880 

0,790 

0.854 

0.645 

0.599 

0.564 

0.510 

0.389 

0.456 

0.316 

0,339 

0.318 

0.273 
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power  spectrums  to  the  original  surface  power  spectrum. 
The  resampling  and  spectra  are  computed  over  the  regions 
for  which  the  resampling  with  thrown-away  rows  and  col¬ 
umns  were  performed  as  described  in  Section  6.2  above. 
Thus  the  resulting  grids  are  54  x  54  and  38  x  38  and  the 
information-containing  portions  of  the  (even)  power  spec¬ 
tra  out  to  the  foldover  frequencies  are  of  lengths  28  and 
20.  When  studying  these  spectra,  a  useful  and  simple 
rule  is  that  a  given  row  N  corresponds  to  a  spatial  fre¬ 
quency  of  (N-l)/M  cycles  per  row  or  column  where  M  is 
either  54  or  38.  Thus,  for  example,  the  numbers  28  and 
20  correspond  respectively  to  (28-l)/54  and  (20-l)/38 
cycles  per  row  or  column  or  exactly  the  foldover  fre¬ 
quency  ratio  of  1  cycle  to  2  samples. 

Both  algorithms  showed  increased  energy  at 
the  lower  frequencies.  The  Jancaitis  algorithm  seems  to 
do  a  little  better  at  retaining  attenuated  frequencies, 
those  for  which  there  is  less  power  in  resampled  grids 
than  in  the  raw  data  grids.  This  advantage  is  always 
less  than  10  percent  and  does  not  translate  itself  into 
noticeably  better  performance  with  regard  to  visual  exam¬ 
ination  of  the  resampled  grids.  Interestingly  enough,  in 
three  cases  the  Akima  algorithm  does  significantly  better 
at  the  very  highest  frequencies  around  the  foldover 
point. 
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6.6 


FINAL  CONCLUSIONS  AND  RECOMMENDATIONS  OF  RESAM- 


PLING  STUDY.  The  two  algorithms  give  roughly  the  same 
results  in  performing  resampling  over  the  data  sets 
provided.  Measured  by  visual  results  and  by  overall 
statistics  there  was  no  way  to  pick  a  best  algorithm.  In 
terms  of  frequency  characteristics  the  Jancaitis 
algorithm  seemed  to  perform  quantitatively  better,  but 
this  did  not  appear  to  have  noticeable  visual  impact. 
Possibly  in  other  data  with  a  higher  frequency  content 
this  would  have  been  more  apparent.  In  dealing  with 
overshoot  the  Akima  algorithm  is  much  superior  in  certain 
situations.  In  terms  of  computational  requirement  the 
Jancaitis  algorithm  performed  significantly  faster  as 
implemented.  . 

Thus  it  would  seem  that  the  question  is 
one  of  trading  the  superior  speed  of  the  Jancaitis 
algorithm  against  the  transition  properties  of  the  Akima 
algorithm.  Unless  it  is  possible  by  skillful  programming 
or  data  handling  to  overcome  the  2.5  to  1  advantage  of 
the  Jancaitis  algorithm,  or  unless  the  overshoot  problem 
is  too  serious  and  too  difficult  to  deal  with  in  other 
ways,  it  would  seem  that  the  Jancaitis  algorithm  provides 
the  better  choice. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


7.1  CONTOUR-TO-GRID.  The  Contour- to-Orid  algorithm 
( CTOG )  has  been  implemented  for  use  at  ETL  and  DMA  in  a 
research  and  testing  environment.  Evaluating  the  grids 
produced  by  CTOG  is  difficult.  Examining  individual  grid 
values  is  time-consuming  and  reveals  only  gross  errors. 
Contour  lines  generated  from  the  grid  can  easily  be  com¬ 
pared  to  contour  line  input,  but  this  introduces  an  un¬ 
certainty  of  the  source  of  an  error  and  provides  little 
information  about  the  value  of  the  grid  elsewhere.  How¬ 
ever,  given  these  constraints,  the  results  from  the  tests 
using  small  subsets  of  the  supplied  data  indicate  that 
the  output  grids  correspond  well  to  the  input  contour 
strings,  ridge  and  drain  lines,  and  lake  boundaries. 

A  few  areas  in  the  test  case  contours  show 
broken  contours,  extraneous  closures,  and  small  bumps. 
These  problems  are  associated  with  small  flat  regions 
where  all  grid  elevations  are  exactly  equal  to  a  contour 
level.  These  result  from  the  coarseness  of  the  digitiz¬ 
ing  increment  relative  to  feature  size  in  the  input 
data.  It  is  recommended  that  scan  step  size  be  smaller 
than  one-half  the  size  of  the  smallest  feature  that  is 
expected  to  be  reconstructed  and  the  CTOG  gridding  incle¬ 
ment  be  approximately  equal  to  this  smallest  feature 
size.  This  guideline  will  result  in  a  very  large  number 
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of  points  to  represent  the  contours,  which  suggests  that 
controlled  thinning  and  perhaps  smoothing  the  input  data 
would  be  advised.  This  step  would  increase  processing 
speed  dramatically  with  little  or  no  loss  of  accuracy. 

7.2  GRID- TO- CONTOUR.  The  Grid- to-Contour  algorithm 
(GTDC)  has  also  been  implemented  for  use  at  ETL  and  DMA. 
The  validity  of  this  algorithm  is  only  slightly  more  evi¬ 
dent  than  for  CTOG.  Although  some  verification  of  con¬ 
tours  from  visual  inspection  of  grid  values  can  be  accom¬ 
plished,  subtle  interpolation  eccentricities  cannot 
easily  be  detected.  Comparison  of  contours  input  to  CTOG 
v tth  those  output  from  GTOC  can  provide  some  assurance, 
and  in  most  parts  of  the  test  areas  the  agreement  is 
high.  The  contour  algorithm  has  some  difficulty  with 
flat  areas  of  height  equal  to  the  contour  level.  The 
algorithm  currently  stops  drawing  when  the  contour  enters 
a  cell  for  which  the  contour  path  is  undefined.  By  digi¬ 
tizing  at  an  appropriate  step  size  for  the  features 
represented,  many  of  the  problem  areas  will  be  elimi¬ 
nated.  For  those  remaining,  a  possible  solution 
would  be  to  choose  a  standard  resolution  for  the  ambigu¬ 
ity,  such  as  always  contouring  slightly  above  a  flat  area 


at  the  contour  level. 


' . 3  SURFACE  SMOOTHING  AND  CONTOUR  GENERALIZATION. 
Two  different  types  of  surface  smoothing  algorithms  have 
been  implemented  and  graphic  results  are  shown  for 
several  variations  of  each  type.  Contour  generalization 
is  provided  by  the  use  of  these  smoothing  techniques  with 
constraints  on  the  amount  and  location  of  smoothing  per¬ 
formed.  The  tests  indicate  that  this  approach  to  contour 
simplification  provides  results  similar  to  direct  line 
simplification  with  some  advantages.  For  areas  with 
dense  contours  the  amount  of  computation  could  be  consid¬ 
erably  less  for  surface  smoothing  than  for  some  types  of 
line  smoothing.  The  process  is  also  more  easily  con¬ 
trolled  to  specifications  normally  used  for  contour  line 
accuracy  and  it  does  not  allow  smoothing  which  violates 
the  meaning  of  the  contour  (allowing  contours  to  cross 
for  example). 

Several  measures  of  filter  convergence  for 
use  with  the  bi-harmonic  filter  have  been  examined.  The 
criterion  of  percent  change  in  surface  curvature  seems  to 
best  represent  the  goal  of  minimizing  curvature  within 
the  constraints  of  the  input  data,  although  the  less 
costly  measure  of  A Z  Max/Range  would  perform  well  for  the 
cases  tested. 
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GRID  RESAMPLING 


7.4  GRID  RESAMPLING.  Several  performance  compari¬ 
sons  between  two  interpolation  algorithms  used  for,  grid 
resampling  have  been  made.  The  two  algorithms  are  the 
Akima  local  bicubic  and  the  parabolic  Jancaitis.  Statis¬ 
tical  comparison  of  interpolation  errors  and  visual  exam¬ 
ination  of  the  output  created  by  the  two  algorithms  show 
that  they  provide  very  similar  performance. 

The  Jancaitis  algorithm  has  a  significant 
speed  advantage  for  the  implementation  scheme  utilized  in 
these  tests.  Other  approaches  which  took  advantage  of 
favorable  relationships  between  input  and  output  grid 
spacing  would  reduce  this  advantage. 

The  Akima  method  has  a  better  ability  to 


model  severe  transitions  without  overshoot.  This  advan¬ 


tage  did  not  show  up  in  the  test  cases  but  is  easily 
quantified.  The  frequency  response  of  the  two  algorithms 
is  comparable,  with  a  slight  advantage  for  Jancaitis. 
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LINEAR  AND  HIGHER  ORDER  FITTING  IN  CONTOUR-TO-GRID 

INTERPOLATION 

In  the  transformation  from  contour  strings  to 
digital  terrain  model  it  is  necessary  to  perform  interpo¬ 
lations  along  elevation  profiles  to  obtain  elevation  es¬ 
timates  at  grid  nodes.  This  is  done  in  two  different 
ways.  One  is  a  simple  linear  scheme  based  on  two  eleva¬ 
tions  of  the  profile  between  which  the  grid  node  is  posi¬ 
tioned.  The  other  involves  the  weighted  average  of  two 
overlapping  quadratic  functions.  This  latter  method  is  a 
natural  adaption  of  the  two  dimensional  Jancaitis  inter¬ 
polation  algorithm  discussed  in  Appendix  F. 

In  order  to  understand  these  two  methods  assume 
that  (Xj ,Ei ) , (X2 »E2 ) » (X3 # E3 )  and  (X4,E4)  are  four  adja¬ 
cent  terrain  profile  pairs  giving  the  relative  position 
along  the  profile  in  the  X  coordinate  and  the  elevation 
in  the  E  coordinate.  Furthermore  assume  that  a  grid  node 
is  located  between  the  second  and  third  elevation  at  the 
relative  position 

p  =  ( X-X2 )/( X3-X2 ) 

For  the  linear  interpolation  method  the  elevation  at  the 
grid  node  will  be  defined  by  the  weighted  average  of  E2 
and  E3  given  by 


Z  =  ( 1-p) E2  +  PE3 


For  the  quadratic  interpolation  scheme  two 
quadratic  functions  will  be  fit  to  the  two  sets  of  data 
pairs  (Xi,Ei),(X2»E2)#(X3,E3)  and  ( X2 , E2 ) , ( X3 , E3  ) ,  and 
(X4,E4).  The  coefficients  for  the  first  quadratic  will 
be  and  Ci  which  solve  the  set  of  linear  equations 
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In  a  similar  manner  the  coefficients  for  the  second 
quadratic  will  be  A2,B2,  and  C2  which  solve  the  set  of 
linear  equations 
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Then  two  estimates  of  the  grid  node  elevation  will  be 
provided  by 

Z1  =  A^2  +  Bjx  +  Ci 
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and 


Z2  =  A2X^  +  B2X  +  C2 

Finally  these  two  estimates  are  averaged  by  use  of  the 
weighting  function 

W(p)  =  (1-pMl-p)  (2p+l) 
to  create  the  combined  estimate 

Z  =  W(p)Zi  +  ( l-W(p) )Z2 
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CONTOURING  ALGORITHM 


This  appendix  contains  a  more  detailed  discus¬ 
sion  about  the  logic  used  in  the  Curve  Priority  Contour- 
ing  program  developed  for  DMA.  The  subject  was  original¬ 
ly  introduced  in  Section  4.0. 

B.l  Contour  Tracing 

All  contours  at  level  C  of  Z(x,y)  are  generated 
by  systematically  examining  each  grid  cell  for  contour 
crossings.  If  a  crossing  is  encountered,  then  the  ini¬ 
tial  point  (x0,y0)  is  computed.  If  (x0,y0)  is  a  point  on 
a  previously  generated  curve,  then  that  curve  is  ignored 
and  another  crossing  is  sought.  When  a  new  curve  is  lo¬ 
cated,  points  along  the  curve  are  generated  for  graphic 
connection  and  display.  When  a  new  point  along  the  curve 
coincides  with  the  curve's  initial  point  (x0,y0)  or  is  on 
the  perimeter  of  the  gridded  area,  then  the  process  tries 
to  locate  another  curve  at  the  same  level  C.  After  all 
grid  cells  have  been  examined,  another  contour  value  is 
selected  and  the  process  repeats.  It  terminates  after 
all  contours  of  each  level  have  been  generated. 

The  algorithm  used  to  detect  a  contour  crossing 
is  very  simple.  If  adjacent  grid  values  on  the  top  or 
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left  aide  of  a  grid  cell  bound  the  contour  level  C,  then 
the  curve  enters  the  cell.  Otherwise,  it  does  not. 
Although  this  may  not  detect  all  types  of  multiple  cross¬ 
ings,  that  is  not  a  problem  since  it  is  most  improbable 
that  this  scheme  will  fail  to  locate  each  curve  somewhere 
within  the  gridded  area. 

The  initial  point  (x0,y0)  for  any  contour  is  on 
the  grid  line  segment  joining  the  two  grid  values  which 
bound  the  contour  value  C.  Points  along  the  curve  con¬ 
sist  of  all  intersections  of  the  curve  with  the  grid 
lines  plus  additional  intermediate  points  as  needed  to 
insure  a  smooth  graphic  representation  when  the  points 
are  connected.  The  intermediate  points  along  a  curve 
consists  of  some  or  all  curve  intersections  with  interme¬ 
diate  grid  lines  laid  over  the  grid  cell  areas  through 
which  the  contour  passes.  An  example  of  an  intermediate 
grid  is  illustrated  in  Figure  B.l.  This  shows  the  two 
types  of  intersection  points  that  are  computed. 

Intermediate  grids  can  be  selected  as  2x2,  5x5, 
9x9,  or  17x17  depending  on  the  accuracy  required  of  the 
final  product.  Intermediate  grid  lines  are  constructed 
so  that  the  exterior  lines  of  the  lattice  coincide  with 
the  original  grid  lines  which  define  the  cell.  Interior 
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Figure  B.L  Two  ‘lypos  ot  (JoriLom  Intel  i>c;ct  ion  Points 
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lines  are  equally  spaced.  The  2x2  cases  uses  a  simple- 
method  in  which  the  interpolation  is  linear  over  the  fout 
triangular  regions  created  by  the  cell  center  and  each  of 
the  four  sides.  The  cell  center  is  assigned  an  elevation 
equal  to  the  verage  of  the  four  nodes.  This  approach 
permits  a  unique  determination  of  the  contour  exit  from  a 
cell.  On  the  other  hand,  the  5x5  grid  essentially  re¬ 
duces  the  grid  spacing  to  1/4  original  spacing  while  com¬ 
puting  intersections  with  a  more  sophisticated  interpola¬ 
tion  scheme.  Similarly  the  9x9  and  17x17  grids  reduce 
the  grid  spacing  by  factors  of  1/8  and  1/16  respectively. 

Curves  are  "followed"  from  one  cell  to  the 
next.  Starting  with  the  initial  point  (x0,y0),  interme¬ 
diate  points  through  the  first  cell  are  computed.  This 
leads  to  an  adjacent  cell  where  another  set  of  intermedi¬ 
ate  points  are  computed.  This  is  continued  until  the 
curve  terminates.  The  entire  sequence  of  points  is  out¬ 
put  for  graphic  connection  and  annotation. 

Intermediate  contour  intersections  are  computed 
by  two  methods.  The  first  is  called  the  Iterative  Method 
and  the  second  is  the  Stepping  Method.  The  Iterative 
Method  is  designed  to  compute  only  those  points  necessary 
to  yield  a  smooth  curve  through  the  cell.  To  describe 
this  process,  momentarily  assume  that  the  curve's  entry 


point  (x0/y0)  and  exit  point  (xe,ye)  are  known.  The  ob¬ 
jective  is  to  compute  a  minimum  number  of  intermediate 


points  ( x^yA ) , ( X2 » y2 ) • • •  through  the  cell  and  maintain  an 
accurate  fit  to  the  true  contour. 

Figure  B.2  illustrates  the  process.  The  dashed 
curve  shows  the  path  of  the  unknown  contour  through  the 


cell  • 

Initially  a  straight 

line  joins 

u0 

*y0> 

and 

( *0  f y e ) 

as  a  first  estimate  of 

the  curve. 

To 

test 

the 

accuracy,  the  displacements  |x0-xe|  and  |y0-ye|  are 
computed.  The  larger,  say  |x0-xe|  is  halved  and  rounded 
to  the  x  coordinate  of  the  closest  intermediate  vertical 
grid  line.  For  example,  with  j  xQ-xe|=  1  and  a  9x9  local 
grid,  the  rounded  x-coordinate  is  Xf=.5.  The  objective 
is  to  compute  y  which  solves 

Z(xf,y)-C  =  0  (B.l) 

where  Xf  is  fixed  and  C  is  the  contour  value.  This  is 
done  iteratively.  An  initial  guess  for  y  is  computed  as 
the  intersection  of  the  straight  line  joining  (x0,y0)  to 
(xe,ye)  and  the  vertical  intermediate  grid  line  x=xf. 
This  yields  yf  which  in  turn  identifies  the  two  interme¬ 
diate  horizontal  grid  lines  which  cross  above  and  below 
yf  at  location  ya  and  y^  respectively  on  x=Xf.  Then 
’ '  x  f  *  y  a )  an,3  z(xf*Yb)  are  computed.  If  these  two  values 
'•M*nd  the  contour  value  C,  then  the  intersection  y*  is 
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(Xc.,Ye) 


Fust  Estimate 

Bisector 


Step  1.  The  first  estimate  is  bisected. 

The  curve  inti  ; sects  the  bisec¬ 
tor  at  tiie  new  point.  The  dis¬ 
placement  tium  the  first  esti¬ 
mate  to  tiie  cu  ve  is  shown. 


Step  2.  The  second  estimate  line  segments 

are  bisected  and  curve  intersections 
computed  to  obtain  two  new  points. 
Tiie  curve  may  now  be  drawn  with 
four  line  segments. 


Figure  B.2  The  Iterative  Metlvod 
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computed  by  inverse  linear  interpolation  between  ya  and 
yb  on  x-xf. 

If  the  values  do  not  bound  the  contour  value, 
then  the  slope  computed  using  Z(Xf,ya)  and  Z(Xf,yb)  is 

used  to  predict  which  direction  to  step  along  x=xf  in 

search  of  the  intersection.  If  the  slope  says  to  move  up 
the  line,  then  yb  is  replaced  by  ya  from  the  previous 

step  and  a  new  ya  is  computed.  This  repeats  until  the 

final  point  (xf«Yf)  on  the  curve  is  calculated.  It  is 
added  to  the  string  of  points  as  (*i»Yi)  to  obtain 

( x0,y0) ' ( X1 'Y1 > >  ( xe,ye >  where  e  =  2. 

If  the  distance  between  the  first  and  last  es¬ 
timates  for  yf  is  less  than  the  spacing  between  local 
grid  lines,  then  no  more  intermediates  are  calculated. 
Otherwise,  the  largest  of  the  four  distance J Xq-x^ | , 

| Yo-yi | » | xl-x2 ) '  and  jyi '-Y2 | is  selected,  halved,  and 
rounded  to  the  coordinate  of  the  nearest  local  grid 
line.  Using  the  procedure  described  above,  this  will  in¬ 
sert  a  new  point  between  (x0,y0)  and  (x^,yj)  or  (xj_,yj_) 
and  ( x2  »Y2 ) ’ 

When  a  new  point  is  inserted  between  two  former 
points  and  when  the  displacement  from  the  initial  to 
final  values  for  the  new  point  is  large  relative  to  the 
local  grid  interval,  new  points  are  necessarily  calcu¬ 
lated  for  the  two  intervals  on  either  side  of  the  new 
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point.  That  is,  if  (xj^y*)  is  a  new  point,  and  the  dis¬ 
placement  is  large,  then  the  intervals  between 
<xk-l'Yk-l>  and  (xk'Yk)  and  between  (xk»Yk)  and 
(xk+l»Yk+l)  are  also  divided  by  new  points.  If  the  new 
point  (Xj,yj)  between  (x^-^.y^-i)  and  (xj^y^)  has  a  small 
displacement,  then  no  more  points  are  added  to  that 
interval.  It  remains  as  ( x^-i  ,yjc-l> » ( x  j  ,y j  ) ,  (  xj^y^) . 
Alternately,  if  the  new  point  (xm,yn)  between  (x^y^)  and 
(xk+l'Yk+l)  has  a  large  displacement,  then  two  more 
points  must  be  computed.  This  would  yield  a  sequence 

( xk'Yk ) *  t xm-l 'Ym-1 )  *  ( xm 'Ym^ •  t xm+l • Ym+1 ^ •  ( xk+l ' Yk+1 ) • 

Note  that  this  process  terminates  with  3  to  5  points 
across  a  grid  cell  if  the  curve  is  nearly  straight.  More 
points  are  generated  when  there  is  considerable  surface 
curvature  across  the  grid  cell. 

For  the  Iterative  Tracing  Method  to  work,  local 
curvature  over  the  area  defining  Z(x,y)  must  be  relative¬ 
ly  small.  Otherwise,  it  is  possible  for  one  of  the  se¬ 
lected  local  grid  lines  to  intersect  the  curve  two  or 
more  times.  This  could  cause  the  algorithm  to  fail.  For 
example,  suppose  x=xf  is  selected  and 
Z(xf,y)-C  =  0 

is  to  be  solved.  If  there  are  two  y  values  which  satisfy 
this  equation,  it  is  possible  that  the  logic  will 


B-8 


“bounce"  back  and  forth  and  not  find  either  one.  Alter¬ 
nately,  it  could  find  the  wrong  one  for  the  current 
situation.  When  the  Iterative  Method  fails,  the 
Stepping  Method  is  used  to  compute  the  intermediate 
points . 

The  Stepping  Method  is  a  micro-level  tracing 
method.  That  is,  suppose  that  a  local  grid,  either  5x5, 
9x9,  or  17x17  is  layed  over  the  cell  through  which  the 
contour  is  known  to  enter  at  (x0,y0).  If  Z(x,y)  is  known 
at  each  node  of  the  local  grid,  following  the  curve 
through  the  local  grid  is  fairly  simple.  Initially,  the 
curve  is  on  the  edge  of  a  local  grid  cell.  It  must  cross 
one  of  the  remaining  3.  Which  one  is  determined  by  find¬ 
ing  the  Z  values  which  bound  C,  the  contour  level.  When 
the  pair  are  found,  the  contour  point  is  computed  by  in¬ 
verse  linear  interpolation.  The  new  point  is  shared  by 
an  adjacent  local  grid  cell.  Therefore,  the  process  of 
checking  the  other  three  sides  is  repeated.  The  stepping 
from  one  local  grid  cell  to  the  next  terminates  when  the 
exit  point  from  a  cell  is  on  one  of  the  original  grid 
lines.  Of  course,  this  leads  to  a  new  grid  cell  and  the 
Iterative  Method  is  evoked  to  try  for  the  next  step  of 
intermediate  points.  The  Stepping  Method  is  used  only 
when  the  Iterative  Method  fails. 
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In  actual  operation,  the  Stepping  Method  does 
not  require  all  values  of  Z(x,y)  on  the  local  grid.  Only 
values  along  the  path  of  the  curve  are  computed  as  they 
are  required.  However,  for  a  17x17  local  grid,  at  least 
15  points  are  computed  if  the  curve  passes  between  oppo¬ 
site  sides  of  the  grid  cell. 

B . 2  Definition  of  the  Local  Elevation  Model 

The  above  approach  requires  a  mathematical  rep¬ 
resentation  E(x,y)  of  the  elevation  over  the  grid  cell  in 
question.  Moreover,  that  equation  must  be  simple  to 
solve  in  both  the  formulation  of  Equation  B.l  and  the 
alternate  formulation 

Z(x,yf)-C  =  0  ( B. 2 ) 
in  which  y  is  fixed  and  x  is  allowed  to  vary.  In  the 
approach  used  here,  the  elevation  is  locally  approximated 
using  16  adjacent  DTM  grid  values. 

Figure  B.3  illustrates  how  these  values  are 
selected.  The  cell  being  processed  is  at  the  center  of 
the  4x4  sub-grid  extracted  from  the  DTM.  This  collection 
of  16  values  can  also  be  partitioned  into  four  3x3  sub¬ 
grids,  all  of  which  have  the  center  cell  in  common.  As 
described  in  Appendix  F,  which  is  patterned  after  materi¬ 
al  by  Jancaitis,  a  quadratic  fit  can  be  made  to  each  of 
the  3x3  sub-grids  to  produce  4  separate  approximations  to 
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the  elevation  over  the  common  center  grid  cell.  These 
are  called  E^(x,y)  through  E4<x,y).  To  produce  a  compos¬ 
ite  E(x,y)  the  four  approximations  are  weighted  and 
averaged  using  the  formula 

4  4 

E( x,y)  =  Wi(x,y)Ei(x,y)/ 

i=l  i=l 

In  the  event  that  one  of  the  quadratic  fits,  E^(x,y), 
cannot  be  constructed,  then  the  corresponding  weight, 
wt(x,y) ,  is  identically  zero.  This  means  that  E(x,y)  is 
usually  the  average  of  four  estimates;  however,  in  spec¬ 
ial  cases  such  as  along  the  edge  or  at  the  corner  of  the 
DTM,  it  may  be  composed  of  fewer  terms. 

Because  the  weighting  formula  is  third  order, 
E(x,y)  is  actually  a  fifth  order  polynomial  in  both  x  and 
y.  Since  solving  Equation  B.l  or  B.2  directly  could  be 
somewhat  time-consuming,  a  further  approximation  is 
made. 

Recall  that  to  define  Equations  B.l  and  B.2  a 
sub-grid  lattice  was  imposed  over  the  cell  being  process¬ 
ed.  If  this  sub-grid  consists  of  17  rows  and  17  columns, 
where  the  exterior  rows  and  columns  coincide  with  the 
four  edges  of  the  DTM  cell,  then  there  are  289  locations 
within  the  cell  where  values  of  E(x,y)  might  be  tabu¬ 
lated.  From  the  tabulated  values  it  is  possible  to  use 


Wj.(  x,y) 
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inverse  linear  interpolation  to  solve  Equation  B.l  or 
B.2.  For  example,  if  the  tabulation  values  along  a  hori¬ 
zontal  sub-grid  line  y=yk  are  labeled  Elx^y^), 


e( xi ,yk)  / 

. . . , e ( x^7 t yk ) , 

and 

if 

two 

adjacent  values 

(E(*i»yk) 

and  E(xi+1,yk) 

,  in 

this 

set 

of  17  elevations 

straddle 

the  contour  L, 

then 

the 

x-coordinate  where  the 

curve  crosses  the  grid  line  y=yk  can  be  computed  using 


x 


xi 


+ 


E(xi,yk)-L 


E(Xi,yk)-E(xi+iyk) 


( xi+l-xi ) 


Obviously,  this  formulation  is  much  simplier  than  solving 
a  5th  order  polynomial  for  its  roots  along  the  line  seg¬ 
ment  y=yk. 

A  similar  formulation  can  be  produced  for  any 
vertical  sub-grid  line  x=Xj.  The  node  values  here  are 
E(xj  »y0) *  E<Xj,yi>,  E(xj,y2),  etc.  If  the  curve  crosses 
between  two  nodes  with  elevation  E(Xj,yk)  and  E(Xj,yk+i) 
on  the  line  x=xj,  then  the  y-coordinate  of  the  intersec¬ 
tion  (Xj,y)  is  computed  by 


y 


yk  + 


E ( Xj , y  k ) — L 
( e ( x j ,yk)-E( x j ,yk+i ) 


(Yk+l-Yk*  • 
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B.3 


INPUT  DATA 


The  input  data  for  CONTRC  is  a  gridded  terrain 
elevation  model  having  NROWS  horizontal  grid  lines  and  N 
COLUMNS  vertical  grid  lines.  The  data  must  be  stored  on 
a  disc  file  columnwise  with  NROWS  entries  per  record. 
The  format  is  sequential  and  binary. 
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appendix  c 

DERIVATION  OF  THE  HI -HARMONIC  OPERATOR 


The  followin«j  derivation  is  along  the  lines  presented 

by  Briggs  17]  and  Young  [8]. 

The  final  result  of  a  yridding  procedure  should  be  a 
grid  model  which  accurately  fits  the  data  and  which  varies 
smoothly  between  the  data  points.  In  mathematical  terms,  this 
smoothness  requirement  may  be  rephrased  as  imposing  a  minimum 
curvature  constraint  upon  the  gridded  array.  This  constraint 
provides  a  way  to  adjust  the  grid  values  to  provide  optimal 
smoothness . 

Let  (x.,y.,z.  . )  be  a  gridded  surface  array.  The 

i  1  i  /  3 

discrete  total  squared  curvature  may  be  defined  as 


where  the  curvature  of  z(x,y)  at  <i,j)  is  defined  as: 


for 

or, 

for 

or, 

for 


c.  =  z .  ,  ,+z.  ,  ,  +  z.  ,+z.  . ,  ,  -  4  z  . 

i,j  1+1,3  1-1,3  i/l-l  i/l+l  i/l 


2  s  i  ^  1-1  and  2  £  j  <  J-l  ; 


c.  .  =  z  .  ,  .  .  +  z  .  .  .  -  2  z  . 

1,3  1+1,3  i"l/l  i/l 


2  <  i  <  1-1  and  j=l  or  j=J 


c.  .  =  z .  .,+z.  .  ,  -  2z .  , 

1,3  1,3+1  i»l-l  i/l 


2  <  j  <,  J-l  and  i=l  or  i  =  I 


C  is  minimized  when 


JC 


r>zi  /  j 


=  0  for  all  i,j 
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Six  cases  must  be  considered: 


1. 

2<i< 1-1 

t 

2-  3- 

1-1 

2. 

2< i< 1-1 

t 

j  =  2 

/ 

3. 

2< i< 1-1 

t 

3-1 

» 

4. 

i  =  2 

t 

3  =  1 

/ 

5. 

i=2 

i 

3  =  2 

t 

6. 

i=  1 

$ 

3  =  1 

. 

All  other  cases  may  be  obtained  by  rotations  of  these  six. 


Consider  first  Case  1.  The  minimization  requirement 
may  be  written  in  the  expanded  form 


<»C 


3  z  . 


a<ci-*-l,j) 


i /  3 


Dz  . 


1 ' 3 


3  Z  . 


1/1 


3  (C .  .  .  ) 

,  1  tL- 

JZi  /  j 


a  (c. 


-izA. 


3zi,i 


3z 


i-3 


Since  z.  .  is  only  contained  in  the  curvature  terms  given,  these 

1  •  J 

five  terms  may  be  expanded  further: 


3(c.Al  .) 

H-l/3 _  =  _ 

3  z  •  .  3  z 


i/l 


— —  J"  z  .  ,  .  +  z  . 

i,j  L  1+2'3 

2  ^  Zi  +  2  ,  j  +  zi,  j 


+  Zi+l,j*-'  +  Zi+1 ,  j-1  “  42i+l 


•J 


J  +  zi+l,j+l+  zi+l,j-l  ”  4zi+l 


■0 


3 (c .  ,  ,) 

3z .  • 

i/l 


=  2 


i.i  [2i'j  *  + 

*  -wj  -  V.,..  *  *i-i.l-i  -  ’ 
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3  (c 


3  z 


i ,  3 


3  z 


=  2 


3  r* _  .  .  ,  +  z.  +  z.  .  -  4z. 

-  iH,]  +  l  t-1,3+1  i  ,]  +  2  1,3  l/3  +  l 

i*3  L  J 

[z  .  ,  ,  .  ,  +  z .  ,  .  +  z  .  ,  +  z  .  .-4z.  -  ,  |  , 

1+1/3+1  l-l, 3*1  1,3+1  1,3  1,3+1 


3 (c .  .  ,  ) 


^ ^  ~  ^ =  -r— ^ fz.  ,  .  +  z  .  ,  +  Z  .+  Z.  .  „  -  4  z  .  ,1 

3zi.j  32i,3  L 1  *• 3-1  ij  io_2  io_1J 

£zi+l,j-l  +  zi-l,3-l  +  zi,j  +  Zi , j-2  “  4zi,j-lj 


=  2 


and 


3  (C  .  .  ) 

_ L l1- 


32ij 


3  z 


- -  I  z  .  ,  ,+z  ,  .  +  z  .  .  +  z .  .  -  4z .  j 

i.iL1*1'3  1'1-1  10+1  1’’"1  10J 

~8  |^Zi  +  l ,  3  f  Zi-1 ,  j  +  Zi,j  +  1  +  Zi,j~l  ~  4zi,jJ 


Accumulating  terms  with  the  same  subscripts  yields; 


20z.  .  +  z.  _  ,+z  ~  ,+z.  +  z  _ 

1,3  1+2,3  1-2,3  1,3+2  i,3-2 


+  2(z.  ,  ...  +  z  .  .  .  .  +  z .  .  ...  +  2... 

l-l, 3  +  1  l-l, 3-1  i+l, 3  +  1  l  +  l, 3-I) 


-  8(z.  .  •  +  z.  .  .  +  z.  ,+z.  ., . )  ~  0 

1-1,3  i+l,3  i,3-l  i,3+l 


(  C  .  2  ) 
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For  Case  2,  - k-4-** -  =  0 

3  z  . 

1,3 
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This  yields, 


W*i.J  *  Zi*2.j  *  *i-2 ,  j  *  *i,j«  "  +  = 


*  2l2l-l,j.l  *  >1.10.1'  -  42i,j-l 


-  8(2.  ,  .  +  2.±.  .  +  Z  .  ,)=0 

1-1,3  i+l,3  i,3+l 


3 (c  ) 

For  Case  3,  - LU—i. 


3  z  .  . 

i,3 


=  0 


3  (c.  .  . )  z 

- — ^-*-2 =  2  Tz  +z  -2z 

3z.  .  *  j  zi-2,j  Zi,j  ZZi-l, 

i ,  j  L 


and 


fiftti ail,2r,  .2  .  2Z 

3zi,j  L  i+2,j  Zi , j  i  +  l, 


This  yields, 

?Zi , j  +  Zi+2 , j  +  zi , j+2  +  zi-l , j+1  +  zi+l,j+l  + 


-4(z.  ,  .  +  z .  .  , .  +  z  .  ,.)=0 

1-1,3  1,3+1  1+1,3 


Similarly  for  Case  4, 


i+1 , j“l 


0 

0 


6zi,j  +  Zi , j+2  +  Zi+l,j+l  +  zi-l,j+l  +  Zi+2,j 


Case  5 


18zi,j  +  2i  ,  3  +  2  +  Zi  +  2  ,  j  +  zx-l  .  j  +  1  +  Zx  +  1,D-1 


+  2z.  .  .  ,  -  8(z  ..  +  z  ...  )  -  4(z  .  +  z .  .  .  )  =  0 

l+l, j+1  i,3+l  i+l,3  i,j-l  1-1,3 


( C .  6  ) 


and  Case  6, 


2Zi ,  3  +  Zi,j  +  2  +  Zi  +  2 ,  j  "  2(zi,3  +  1  +  Zi  +  1,3>  “  ° 


( C  .  7  ) 


Equations  C.2  through  C.7  define  a  system  of  linear 
equations,  one  for  each  grid  position,  which  relate  each  Z^fj 
to  a  pattern  of  surrounding  grid  values.  In  matrix  form  the 
system  can  be  represented  by 


Wz  =  0 


where  z  is  a  vector  consisting  of  ail  grid  entries  while  W  is  a 
square  matrix  that  is  mostly  zero  filled.  The  main  diagonal  is 
unity  with  six  minor  non-zero  diagonals  on  either  side  of  the 
main  diagonal.  All  other  entries  are  zero. 

It  is  possible  to  solve  this  system  of  equations  using 
classical  Gaussian  elimination.  However,  for  many  applications, 
the  number  of  grid  entries  makes  that  impractical.  Therefore, 
iterative  means  are  used  to  "solve"  C.8. 


There  are  various  iterative  techniques  for  solving 
large  systems  of  linear  equations.  The  method  well-suited  to 
this  specific  problem  is  called  Successive  Overrelaxation  or 
Extrapolated  Leiberman.  The  iteration  equation  is 


Zi,j(n)  =  CO  )wkZk(n-l>  +  )wmZm<n> 

(l-0))zi( 


(  C .  9  ) 
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Where  n  denotes  the  iteration  index,  the  computations  are  applied 
systematically  to  the  grid  values  so  that  each  current  Is 

replaced  by  an  updated  Zj^j. 

The  initial  Zj^j's  denoted  by  2  £  #  j  ^  °  ^ ,  must  be 
established  prior  to  the  iterative  update  process  and  they  should 
be  as  close  to  the  final  solution  as  possible  to  minimize  the 
number  of  iteration. 

The  first  summation  in  C.9  spans  the  set  of  Z^j's 
that  have  not  been  modified  by  the  current  update  while  the  sec¬ 
ond  summation  spans  those  that  have  been  adjusted.  For  example, 
if  the  iteration  proceed  down  the  grid  columns  and  from  left-to- 

right,  then  terms  in  the  first  sum  are  always  below  and  right  of 

the  center  while  terms  in  the  second  sum  are  above  and 

left  of  Zj^j. 

The  (J  term  in  C.9  is  called  the  relaxation  factor. 

WhenOJ  =  l»  the  grid  values  are  "relaxed"  by  the  iterations.  That 

is,  the  center  term  Z^(j  is  simply  replaced  by  a  weighted  com¬ 
bination  of  up  to  12  surrounding  grid  values.  The  specific 
weights  depend  on  the  position  of  Zj#j  in  the  matrix  as  de¬ 
scribed  by  Cases  1  through  6. 

For  (0  >1,  over-relaxation  results.  That  is,  the  new 
Zifj  that  would  be  obtained  with  <J=1  is  linearly  extrapolated 
forward.  For  example,  if  and  Z  ^  r  j  ^ n  ^  are  two 

successive  iterations  produced  with  GJ  =  1  and  both  are  moving  in 
the  direction  of  the  solution,  then  the  rate  of  convergence  can 
be  accelerated  by  extrapolating  forward. 

Theoretically,  there  is  an  optimum  value  for  (jJ  that  is 
somewhere  between  1  and  2.  It  can  be  shown  that  use  of  smaller 
than  the  optimum  value  slightly  reduces  the  convergence  rate. 


C-6 


Also,  use  of  (jJ  greater  than  the  optimum  value  yields  instabili¬ 
ties  and  possible  divergence.  Since  the  theoretical  (jJ  is  diffi¬ 
cult  to  establish,  a  fixed  CJ  is  typically  used.  For  most  situa- 
tions  GJ  =  1  •  3  is  satisfactory  and  values  as  large  as(t)  =  1.5  might 
work  but  frequently  cause  instabilities  in  the  process.  Values 
less  than  1.3  simply  reduce  the  rate  of  convergence  below  that 
observed  by  the  emperically  selectedCJ=l.  3. 

Not  all  grid  values  are  recomputed  by  C.9.  Fixed  or 
known  grid  values  are  skipped  as  the  adjustments  are  applied  sys¬ 
tematically  to  the  matrix.  Grid  entries  between  fixed  values  are 
adjusted  to  achieve  the  minimum  curvature  results.  Fixed  grid 
values  are  analogous  to  boundary  values  in  differential  equation 
problems. 

The  iterations  of  C.9  with  the  grid  are  applied  system¬ 
atically  in  groups  of  four.  The  first  of  four  starts  at  the  up¬ 
per  left  matrix  position  and  works  down  the  columns  and  from  left 
to  right  from  column-to-column.  This  results  in  a  spreading  of 
information  content  in  fixed  values  toward  the  lower  right  corner 
of  the  matrix.  To  counter  this  effect,  the  second  iteration  is 
applied  in  exact  reverse  order  to  the  first.  The  first  and  sec¬ 
ond  iterations  now  distribute  information  in  diagonal  patterns 
within  the  gridded  area.  This  is  countered  by  iteration  three 
and  four.  The  third  starts  at  the  lower  left  corner  and  works  to 
the  upper  right  corner  while  the  fourth  backtracks  on  the  third. 
The  results  are  a  uniform  distribution  of  information  away  from 
any  point  on  the  grid. 
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APPENDIX  D 

LEAST  SQUARES  I NTERPOLAT 1  ON  IN  TWO  DIMENSIONS 


Consider  a  set  of  data  (x^,y^,2^)  to  which  it  is 
desired  to  fit  the  relation  of  a  plane 

Z(x,y)  =  A(x-Xq)  +  B (y-y  )  +  C  .  (D.l 

Assume  that  the  (x,y)  values  are  precise!,  all  uncertainly  is  in 
the  2  values.  Define  the  squared  difference  to  be  minimized 

N  2 

E  =  Z‘°i[2rz(xx'yi)] 

i  =  l 

where  N  is  the  number  of  points  tc  be  interpolated  to  (x  ,y  ) 

o  o 

and  the  are  the  weights  associated  with  each  of  these  points. 

The  equirement  dE=0  is  imposed  to  determine  the 
minimum  E,  wh-ioh  implies 

N 

If  =  °  =  2ZWi[Zi-A(xi"xo)"B(yi-yo)"C](xi-xo) 
i  =  l 

N 

H  “  0  =  2Zaii[Zi"A(xi"Xo)'B(yi‘yo)  'C](yi'yo)  (D-2) 

i=l 

N 

S  =  0  =  2Zwi[zrA(xi'xo)+B(yi'y0)-c]  • 
i*l 
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(D.3) 


AEw .  (x . -x  )  (v  .  -  y  )  +B;  ■  (y  -y  )  +CEi,> .  (y  .  -y  ) 
1  1  o  Ji  o  >  1  1  1  o  t  J 1  1  o 

i  li 


-  }:ta,izi(yi~yo) 
1 


Alu).  (x-x  )+BEui.  (y.-y  ) +C  Eu-.  =  Lu  .  z  . 
lio  l  1 1  o  1  11 

l  l  li 

For  notational  convenience,  define, 

S  =  L(d  .  (x  -x  )  ^ 

XX  1  1  o 

1 


S  =  Eui .  (X  .  -X  )  (y  .  -y  ) 
xy  .  i  i  o  -*10 


S  »  Eu. (x.-x  ) 
x  ^  1  1  o 


s  =  Eu> .  (y .  -y  ) 
yy  o 


S  =  Eo) .  (y  .  -y  ) 
v  i  ■*  i  o 

J  l 


S  =  Eui.  (x.-x  )z. 
XZ  ^  1  1  O  1 


S  =  Ecd  .  (y  -y  )  z  . 
yz  .  i  ■*  i  ■*  o  l 


S  =  Lai .  z  . 

z  .11 

l 


S0  "  Za)i 

l 


Then  the  normal  equations  may  be  written  in  the  3x3  form 


(D .  4 ) 
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& t  a  grid  node  to  be  interpolated#  x~x0  an<^  y-yo*  Since 
Z (xQ,yo) =C ,  only  C  need  be  solved  for.  The  result  is 
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APPENDIX  E 

WEIGHTING  FUNCTIONS  FOR  DIGITAL  TERRAIN  MODELING 


The  weighting  functions  used  for  cartographic 
tasks  should  be  appropriate  for  the  type  of  data  to  be 
interpolated  and  the  use  intended.  Two  weighting  func¬ 
tions  which  are  used  in  contour  to  grid  interpolation  and 
grid  filtering  are  strictly  dependent  on  relative  dist¬ 
ances.  The  first,  which  decreases  rapidly  with  distance, 
is  shown  in  Figure  E.l(a).  It  is  of  the  form 

w(r,R)  =  (1  -r/R) 2/( r/R ) 2 

where  r  is  the  distance  from  the  grid  node  and  R  is 
equivalent  to  the  maximum  search  distance  about  the  grid 
node.  This  weighting  function  is  suitable  for  data  that 
are  essentially  free  of  random  behavior  or  observational 
error.  Such  data  might  be  called  "deterministic". 

A  second  distance-dependent  weighting  function 
which  decreases  slowly  with  distance  is  shown  in  Figure 
E.l(b).  It  is  of  the  form 

w ( r , R)  =  (  l-r/R)2<l+2r/R) 

where  r  and  R  remain  as  defined  in  the  previous  example. 
This  function  is  desirable  when  working  with  data  which 
have  a  random  or  statistical  component  since  the  random 
component  tends  to  be  averaged  out  by  assigning  roughly 
equivalent  weight  to  many  inputs. 

A  third  weighting  function  which  is  of  particu¬ 
lar  use  in  contour  to  grid  interpolation  is  of  the  form 

w(Zl,Z2,r,R)  =  ( ABS( Zl-Z2)/(r/R) ) 
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Figure  E.l.  Two  types  of  weighting  functions 
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where  Z1  and  Z2  are  contour  elevations  on  either  side  of 
grid  node  along  a  profile,  R  is  the  maximum  search  dist¬ 
ance  about  the  grid  node,  and  r  is  the  current  search 
distance.  This  function  weights  estimates  based  on  two 
elevations  of  similar  magnitude  very  small.  Thus  grid 
nodes  which  are  interior  to  drains  and  ridges  and  thus 
surrounded  on  at  least  three  sides  by  a  single  contour 
show  elevation  variation  based  on  search  lines  which 
cross  a  contour  on  the  fourth  "free"  side. 


E-3 


APPENDIX  F 


APPENDIX  F 

THE  JANCAITIS  INTERPOLATION  ALGORITHM 


The  following  is  a  heuristic  development  of 
ZYCOR's  adaption  of  the  Jancaitis  interpolation  algo¬ 
rithm.  An  original  paper  by  Jancai  n  and  Junkins  dis¬ 
cusses  the  derivation  of  the  weighting  functions  used  in 
the  algorithm  and  the  motivation  for  the  blending  of  re¬ 
sults  obtained  for  four  overlapping  sub-grids  (Reference 
2).  The  use  of  parabolic  fits  within  each  3x3  sub-grid 
is  a  concept  developed  by  ZYCOR. 

The  purpose  of  this  algorithm  is  to  interpolate 
elevation  at  arbitrary  (x,y)  locations  from  elevation 
values  represented  in  a  matrix  or  gridded  form.  The 
algorithm  assumes  that  the  (x,y)  location  where  elevation 
is  to  be  interpolated  is  within  the  gridded  area,  the 
grid  lines  are  equally  spaced  (square  grid  cells),  and 
that  all  grid  nodes  have  an  assigned  elevation  value. 

An  interpolated  value  at  (x,y)  is  computed  from 
16  elevation  grid  values.  The  16  values  are  at  the  nodes 
of  a  4x4  sub-grid  where  the  center  cell  area  contains  the 
(x,y)  location.  The  4x4  sub-grid  is  further  decomposed 
into  four  3x3  sub-grids.  Figure  F.l  shows  the  4x4  sub¬ 
grid  and  how  the  four  3x3  sub-grids  are  extracted. 

This  arrangement  of  4x4  nodes  and  3x3  nodes  was 
selected  to  simplify  the  algorithm's  implementation. 
Note  that  the  central  cell  of  the  4x4  containing  (x,y)  is 
always  in  the  lower  right  cell  position  of  each  3x3. 
This  means  that  one  algorithm  can  be  used  to  interpolate 
the  elevation  for  each  of  the  four  3x3  sub-grid  areas. 
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The  center  cell  bounded  by  nodes 
6,  10,  11,  and  7  is  area  over 
which  interpolation  can  be  per¬ 
formed.  Each  of  the  four  3x3  sub¬ 
grids  containing  this  cell  is  ro¬ 
tated  to  place  it  in  the  lower  right 
corner. 


The  3x3  Sub-grid  Node  Ordering 

The  values  selected  for  node  posi¬ 
tion  1  is  the  entry  at  position  1, 
13,  16,  or  4  of  the  4x4  sub-grid 
depending  on  Case=l,2,3,  or  4. 
Similarly  for  the  other  8  node 
positions.  Interpolation  is 
limited  to  (p,q)  in  the  lower 
right  cell  bounded  by  5,8,9, 
and  6 . 


4x4  nodes  selected  for  the 
upper  left  3x3 


Case 

k=2 


(P2'<32) 


4x4  nodes  selected  for  the 
right  3x3 


Case 
k=  3 


P3rq 


3 


) 


4x4  nodes  selected  for  the 
lower  right  3x3 


Case 
k  =  4 


(P4'<34> 


4x4  nodes  selected  for  the 
left  3x3 


Figure  F.l 


Node  Lay-Outs  for  4x4  and  3x3  Sub-Grids 


All  that  is  required  is  a  simple  mapping  of  (x,y)  to  the 
proper  positions  for  each  3x3. 

The  actual  location  (x,y)  within  the  elevations 
grid  is  converted  to  a  relative  location  within  a  unit 
central  cell  of  the  4x4  using  the  relationships 

xR  =  (x-xQ)/  x  and  yR  =  (yQ-y)/  y  (F.l) 
where  (x0,y0)  is  the  actual  coordinate  of  the  upper  left 
corner  of  the  center  cell.  This  is  the  location  of  the 
node  labeled  6  in  Figure  F.l.  Note  that  0<xr<1  and 

OiyRl1* 

The  following  table  gives  the  relative  location 
(p,q)  of  (  xR'Yr)  for  each  of  the  four  3x3  sub-grid 
cases . 


Table  F.l 

Relationship  between  (p^q^)  and  (xR,yR)  for  3x3  case 


3x3 

Case 

Pk= 

9k= 

k=l 

xR=(x-x0)/  X 

yR=(yo-y>/  y 

k=2 

1-Xr 

vr 

k=3 

i-yR 

1-Xr 

k=4 

XR 

i-yR 

If  the  node  positions  of  each  3x3  are 
re-labeled  for  convenience  as  shown  at  the  bottom  of 
Figure  F.l,  then  the  index  relationship  between  the 
elevations  node  position  in  the  4x4  grid  and  the  four  3x3 
grids  is  given  by  the  following  table. 
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Table  F.2 

Relationship  Between  Node  Indices  of  the  General  3x3  and 

the  4x4  Sub-grid 


3x3 

Node 

Indices 

4x4  Node  Indices 

Case  1 

Case  2 

Case  3 

Case  4 

fcr  ^ 

1 

13 

16 

4 

BSllBi 

2 

14 

12 

3 

3 

15 

8 

2 

■ 

5 

9 

15 

8 

5 

6 

11 

7 

6 

7 

11 

7 

6 

7 

9 

5 

14 

12 

8 

6 

10 

11 

11 

7 

6 

10 

Over  each  of  the  3x3  sub-grids  a  bi-quadratic 
surface  fit  is  made  to  interpolate  a  value  at  (p,q)  where 
(p,q)  comes  from  the  relations  in  Table  F.l  The  inter¬ 
polation  at  ( p , q )  is  performed  by  first  interpolating  in 
the  p-direction.  Figure  F.2  shows  where  E^ ,  E2  and  E3 
are  computed.  The  equations  for  computing  each  of  E^,  E2 
and  E3  is  given  by 

E  i  (  p )  =  B.U?"1  J  ej  +  (l-p2)ej  +  3  +  B-lBll  )&j+6  • 

2  2  (F.2) 

where  the  e^'s  are  the  elevations  in  the  3x3  grid.  Now, 
the  3  E £ '  s  are  interpolated  in  the  q-direction  to  get 
E(p,q).  The  equation  is 

E(p,q)  =  ‘lliilli  E£  (  p)  +  ( 1-q  )  2E2  (  p ) 

+  ElEliJ  E  3  (  p )  .  (  F .  3  ) 
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E,  (p , e 


E0  (p,e 


(P.q)  E(p,q)=E(q,E 


Steps  for  computing  E(p,q) 

1. 

Compute 

E^  using  p,  e^. 

e  4  ' 

and 

e7 

2. 

Compute 

E2  using  p,  e2. 

e  5  ' 

and 

e  8 

3. 

Compute 

using  p,  e-j, 

e6  ' 

and 

e9 

4. 

Compute 

E(p,q)  using  q, 

V 

L2' 

and  E^. 

Figure  F. 

2.  Locations  of 
Computing  E(p 

Elevation  Data  and  Steps  for 

,q) 

With  this  equation  and  the  scheme  for  collect¬ 
ing  the  e^'s  for  each  3x3  it  is  possible  to  compute  the 
elevation  at  (p,q)  for  each  of  the  four  3x3  sub-grids. 


Under  most  circumstances  there  is  little  re-.son 
to  expect  the  four  estimates  to  agree  in  value.  This 
situation  is  resolved  by  averaging  the  four  estimates 
using 


4  q 

y)  =  ^(p^qj^E^Pi^q*)/  yw(pk,qk)  (F.4) 

k  =  l  k=*L 


where  (x,y)  is  related  to  the  (pk,qk)  locations  through 
Table  F.l.  The  weight  function  w{p,q)  has  the  general 


form 


w(p,q)  =  (l-p)2(2p+l) (l-q)2(2q+l)  (F.5) 


where  0£p£l  and  0<q<l.  Note  that  w(p,q)  is  third  order 
in  both  p  and  q  while  E(p,q)  is  second  order  in  both  p 
and  q.  Therefore,  E(x,y)  is  of  the  fifth  order  in  both  p 
and  q  over  the  central  grid  cell. 

The  smooth  nature  of  this  interpolating  scheme 
comes  from  the  blending  of  individual  3x3  interpolated 
values  with  the  weight  function  w(p,q).  Note  that  when 
(xr»Yr)  is  at  the  node  position  6  of  the  4x4,  i.e., 

when  (xR,yR)=(0,0) ,  that 

w(Pi,qi)  =  1 

W(P2'<?2)  =  0  *  (F.6) 

w(P3,q3>  =  0  ,  and 

w(P4»q4>  =  o 
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(P.7) 


Hence  E(x,y)  reduces  to 

E(x,y)  -  E1(0,0) 

which  means  that  the  interpolated  value  is  equal  to  the 
value  at  node  6.  Similar  situations  occur  for  the  other 
three  corner  nodes.  When  (xfy)  is  at  the  center  of  the 
central  grid  cell,  all  Pi's*.5  and  all  qi's».5.  Also  at 
the  center  the  weight  function 

w(  .5,  .5)  -  .25.  ( F. 8 ) 

Therefore,  E(x,y)  reduces  to 

4 

E(x,y)  *  1/4  £  Ek(.5,.5)  .  (F.9) 

k=l 

Since  w(x,y<  is  smooth  over  the  central  cell 
area,  the  weighting  scheme  smoothly  blends  the  variations 
in  interpolated  elevations  from  one  3x3  with  adjacent 
3x3’ s.  Consequently,  an  overall  smooth  interpolation 
scheme  is  obtained. 

When  computing  irregularly  spaced  elevation 
values,  it  is  probably  necessary  to  compute  E(p,q)  and 
w(p,q)  for  each  of  the  four  3x3  sub-grid  areas.  However, 
when  the  algorithm  is  used  to  compute  values  at  the  same 
relative  locations  within  many  grid  cells,  the  number  of 
computations  can  be  reduced  by  tabulating  the  weights  and 
coefficients  of  the  Ek  functions. 

This  requires  four  tables.  One  table  has  the 

values  of 

w ( t)  *  (l-t)2(2t+l)  (F.10) 
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II 


at  values  of  t  equal  to  all  values  of  p  and  q  where 
interpolations  are  made.  From  w(t),  the  weights  w(pfq) 
can  be  computed  with  one  multiplication.  The  other  three 
tables  contain  values  of 

t(t-l)/2 

(1-t2)  ,  and  (F.ll) 
t(t+l)/2 

at  the  same  values  of  t  as  used  to  tabulated  w(t).  With 
these  coefficients  the  polynomials  Ei(p»q)»  E2(p»q)» 
etc.  can  be  computed  with  12  multiplications.  Finally, 
the  weighting  requires  4  more  multiplications.  In  all, 
17  multiplications  are  required  to  compute  E(x,y)  when 
the  interpolating  locations  are  always  the  same  within 
central  grid  cells.  This  compares  with  40  multiplica¬ 
tions  if  E(x,y)  is  computed  without  pre-computed  and  tab¬ 
ulated  weights  and  coefficients. 


F-8 


APPENDIX  G 

THE  AKIMA  BIVARIATE  I NT ERP RELATION  ALGORITHM 


The  Akima  algorithm  (Reference  3  and  4)  is  a 
local  interpolating  technique  in  that  elevations  can  be 
interpolated  by  extracting  a  small  area  of  the  grid 
around  the  location  to  be  interpolated.  The  interpolated 
value  depends  only  on  the  elevations  in  the  small  sub¬ 
grid  and  not  on  values  beyond  the  sub-grid.  The  grid 
values  that  are  required  to  interpolate  within  a  specific 
grid  cell  are  illustrated  in  Figure  G.l.  Here  the  center 
cell  of  the  configuration  contains  the  (x,y)  point  where 
interpolation  is  desired.  The  surrounding  grid  values 
are  required  to  construct  the  interpolating  function  over 
the  central  grid  cell. 

For  convenience  the  (x,y)  location  in  the  cen¬ 
tral  grid  cell  is  converted  to  a  relative  coordinate 
( p , q ) .  Then  all  cells  can  be  treated  as  if  their  sides 
have  unit  length  and  0<p_<l  and  (Kq£l . 

The  interpolating  function  over  the  central 
cell  with  this  algorithm  is  a  bi-cubic  polynomial  of  the 
form  3  3 

e(P.q)  -  £  £  Ri.jpiclj  <G-n 

i  =  0  j  =0 

where  the  16  coefficients  a1,0'  etc.  remain  to 
be  established.  This  requires  at  least  16  elevation  val¬ 
ues  from  the  provided  grid.  Note  that  24  values  are  ac¬ 
tually  marked  in  Figure  G.l.  The  additional  8  values  are 
used  in  computing  weights  which  provide  continuity  and 
smoothness  from  central  cell-to-cell . 
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Figure  G.l.  Pattern  of  elevation  grid  values  required 
to  interpolate  within  the  central  cell 
using  the  Akima  Algorithm 

Only  those  numbered  cells  are  used.  There  are  a  total 
of  24  values  employed  for  each  interpolation. 
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The  16  coefficients  of  the  polynomial  are  com¬ 
puted  using  values  of 

e  ,  ex , Gy ,  and  exy 

at  the  four  corners  of  the  central  grid  cell.  Here  e  is 
the  elevation,  while  ex  and  ey  are  the  partial  deriva¬ 
tives  with  respect  to  x  and  y  and  eXy  is  the  cross  par¬ 
tial  derivative.  These  four  quantities  at  each  of  four 
corners  provide  16  knowns  to  establish  the  16  unknown  co¬ 
efficients. 

The  partial  derivatives  ex,  ey,  and  eXy  for 
each  corner  node  of  the  central  cell  are  computed  using  a 
pattern  of  adjacent  elevation  values.  Figure  G.2  illus¬ 
trates  the  grid  positions  that  are  required.  Here  the 
center  grid  position  (3,3)  corresponds  to  one  of  the  four 
corner  grid  positions  for  the  center  cell.  The  computa¬ 
tions  described  are  repeated  for  each  of  the  four  cor¬ 
ners  . 

The  partial  ex(3,3)  is  computed  as  the  weighted 
average  of  the  left  and  right  difference  equations  for 
e( 3 , 3 ) .  That  is 

ex  ( 3  «  3  )  =  (wX/+ex,_+wx,-eXr  +  )/ 

(wx,-+wx,+)  (G.2) 

where 

eX , -  *  ( e3 , 3~e2 , 3 ) /Ax  *  and 

eX,+  =  <e4,3“e3,3)/Ax 

Since  ax*  Ay  =  l  by  assumption,  division  by  Ax  and  Ay  is 
omitted  in  subsequent  difference  equations. 

The  weights  applied  to  the  left  and  right  par¬ 
tial  estimates  are  proportional  to  the  curvatures  esti¬ 
mated  at  nodes  (3,3)  and  (4,3)  respectively.  The  weights 
are  given  by 
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Indexing  used  to  identity  the  13  elevations 
values  required  to  compute  e  x  ,  e  ,  and  e 
for  each  corner  grid  node.  x  Y 


Indexing  of  nodes  used  to  compute 
and  the  paramaters  involved  in  e  . 

A 


Figure  G.2.  Layout  and  indexing  of  eleva¬ 
tion  grid  values  used  to  com¬ 
pute  ex,  ey ,  and  ex,y 


"X,-  *  |el,3'2e2,3+e3,3 |  *  and 

WX ,  +  "  |e3,3"2e4,3+e5,3 |  •  (G-3> 

In  the  special  event  that  both  *xt-  and  WX»+  are  ^°th 
zero,  then  ex(3,3)  is  computed  by 

ex(3,3)  *  1/2 ( ex , -+ex , + ^  (G.4) 

A  similar  set  of  equations  are  used  to  compute 
ey(3,3).  These  are 

eY  ( 3 , 3  )  *  (wy,+eY,-+wY,-eY,+)/ 

(wY#_+wY#+)  (G.5) 

where  «y,-  and  eY/  +  are  the  difference  equations 
approaching  from  the  top  and  bottom  respectively.  The 
weights  are  given  by 

WY, -  *  | e3 , l”2e3 , 2+e3 , 3 |  '  and 

=  |e3,3"2e3,4+e3,5 1  •  (G.6) 
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The  cross  partial  at  (3,3)  is  somewhat  more  com¬ 
plex.  It  is  estimated  by 


ex ,y ^  ^ ' 3 ) 


-<e3,3-e2,3-e3,2+e2,2) 


+ 


wy ,  + '  e3 , 4‘ 


+ 


w 


X,+ 


wy,-(e4,3'e3,3-e4,2+e3f2) 


+ 


Wy  >  +  (64,4"e3,4~64, 3+e3 , 3 ) 

^Wx,-+Wx,+) (wy,-+wy,+)  •  ( G . 7 ) 

Again,  the  special  case  of  wXf_=wXf+=0  and/or 
Wy#_=Wy>+*0  is  treated  by  setting  the  zero  weights  to 
unity. 

Note  that  when  wx t _=wx t +=wy t _=wy f +=1/2 , 

these  rather  formidable  weighted  difference  equations  re¬ 
duce  to 

ex(3,3)  =  1/2 ( 64 f 3-62 t 3 )  ,  (G.8 

ey(3,3)  =  1/2 ( e3  #  4~e3  f  2 )  *  and 

ex,y(3,3)  =  l/2(e2f2~e2,4+e4,4_e4,2^  » 

which  are  widely  recognized  as  standard  difference  equa¬ 
tions  for  estimating  these  three  partial  derivatives. 

Interpolating  based  on  this  scheme  must  be  con¬ 
tinuous  between  adjacent  central  interpolating  cells  since 
the  same  parameters  are  used  to  construct  the  polynomial 
along  the  interface  lines  between  grid  cells.  For  example. 
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Figure  G.3  shows  6  grid  elevation  nodes.  Along  the  hori¬ 
zontal  grid  line  between  nodes  C  and  D  the  elevations  and 
partials  e^ ,  eD'eX,C»  and  eX,D  are  the  same  whether 
they  are  computed  for  the  central  cell  A,B,C,D  or  for  the 
cell  C,D,E,F.  Since  there  is  only  one  cubic  which  fits 
through  D  and  C  with  values  and  derivatives  given  at  D  and 
C,  the  interface  between  the  two  cells  must  be  continuous. 
Similarly,  since  the  cross  partials  are  computed  so  that 
eX/y=ey/X  the  same  reasoning  implies  that  ey  is  con¬ 
tinuous  across  the  interface.  Thus,  interpolated  values 
will  vary  smoothly  as  the  interpolating  locations  (x,y) 
move  from  one  central  grid  cell  across  the  grid  lines  into 
adjacent  cells. 

The  number  of  operations  required  to  interpolate 
a  single  elevation  with  this  method  can  be  quite  signifi¬ 
cant.  An  optimized  version  of  the  algorithm  designed  to 
use  as  input  the  unit  square  grids  which  are  typical  of  DMA 
applications  requires  78  floating  point  multiplies  and  139 
floating  point  additions  for  each  interpolation.  Thus  any 
method  of  implementation  which  can  improve  the  speed  of  the 
algorithm  is  of  interest. 

The  complexity  of  the  algorithm  is  largly  attrib¬ 
utable  to  the  great  number  of  calculations  necessary  to  ob¬ 
tain  estimates  of  the  various  partial  derivatives  required 
at  each  of  the  four  input  grid  nodes  surrounding  the  point 
at  which  it  is  desired  to  carry  out  the  interpolation. 
Given  the  partial  derivatives,  calculation  of  the  16  coef¬ 
ficients  of  the  bi-cubic  polynomial  used  in  the  interpola¬ 
tion  and  evaluation  of  the  polynomial  itself  requires  only 
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22  floating  point  multiplications  and  64  floating  point 
additions.  This  suggests  that  if  the  interpolation  points 
can  be  ordered  relative  to  the  input  grid  in  such  a  manner 
that  computed  partial  derivatives  can  be  saved  and  reused, 
a  considerable  savings  in  computation  time  could  be 
realized. 

Take,  for  example,  a  situation  where  resampling 
is  at  roughly  the  same  density  as  the  input  grid  and  the 
new  grid  points  have  been  sorted  and  ordered  nicely  so  that 
the  old  grid  can  be  stepped  through  a  column  at  a  time. 
For  each  cell  in  a  particulr  column  partial  derivatives  for 
one  side  will  already  have  been  calculated  when  the  resam¬ 
pling  was  performed  for  the  previous  column.  Also,  both 
sets  of  partials  for  the  top  of  the  cell  would  have  been 
required  to  interpolate  over  the  previous  cell  in  the  same 
column.  Thus,  in  an  average  situation,  only  one  set  of 
partial  derivatives  for  a  single  grid  node  would  need  to  be 
calculated  for  each  interpolation.  Taking  into  account  the 
numbers  mentioned  above  for  floating  point  operations , th is 
suggests  the  possibility  of  reducing  the  computation  time 
for  the  interpolation  operation  by  more  than  50%. 
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